The direct multi-region boundary element method (BEM) is developed to simulate the interaction of ultrasonic waves with submerged finite elastic bodies. This technique is further specialized to model fluid-solid half-space problems with fluid-filled indentations and air bubbles on the solid surface. These problems usually occur in ultrasonic immersion testing and are of particular interest to the NDE community. The problems are formulated in the Fourier (frequency) domain and are described by the three dimensional acoustic and elastodynamic boundary integral equations (BIE) with pressure and displacement serving as primary variables. The techniques developed are general and may be used with any kind of incident wave, although plane waves are used in all numerical experiments. The equations governing the acoustic region are first converted mathematically to equations like those of an elastodynamic region and the resulting two elastodynamic regions are coupled and solved for the interface displacements. Tractions, pressures, and pressure gradients are subsequently computed using these displacements and the interface conditions. Numerical results are obtained for some typical problems and justified using reciprocity relations and by comparison with solutions available for the half-space elastodynamic problem.