This paper is concerned with the three dimensional interaction and coalescence of two bubbles subject to buoyancy and the dynamics of the subsequent joined bubble using the boundary integral method (BIM). An improved density potential method is implemented to control the mesh quality. It helps to avoid the numerical instabilities which occur after coalescence. Numerical convergence tests are conducted in terms of mesh sizes and time steps. The 3D numerical model agrees well with an axisymmetric BIM model for axisymmetric cases as well as experimental results captured by high-speed camera. The bubble jetting, interaction and coalescence of the two bubbles depend on the maximum bubble radii, the centre distance between two bubbles at inception and the angle β between the centre line and the direction of buoyancy. We investigate coalescence of two bubbles for β = 0, π/4, and π/2, respectively, and at various centre distances at inception. Numerical results presented include the bubble and jet shapes, the velocity and pressure fields surrounding the bubbles, as well as the time histories of bubble volumes, jet velocities and positions of centroid of the bubble system.