Why a Shout Above Water Never Reaches Below — Acoustics Across a Diffuse Interface

Call a friend from the edge of a pool and the swimmer underwater barely hears it. Roughly 99.9% of airborne sound bounces off the surface. The same split drives ocean acoustics, medical ultrasound, and underwater noise. Following Ballout et al. (2025), these notes redo reflection, transmission, and refraction with a 1D acoustic solver.

Raise the speed and density ratios and $z_2/z_1$ opens up. The reflected pulse at the center interface grows, and the transmitted pulse shrinks.

Match the two impedances 1:1 and the pulse crosses. The reflection vanishes. An air–water jump leaves only the reflected pulse on the left.

Abbas Ballout et al., "Acoustic Propagation/Refraction Through Diffuse Interface Models", arXiv:2504.01727 (2025). They adjust the weak-compressibility term of an incompressible Navier–Stokes/Cahn–Hilliard system so two phases carry different sound speeds, and the wave rides those speeds.

Acoustic Impedance — How a Medium Resists Sound#

The resistance a medium offers to a passing sound wave is its acoustic impedance (density times sound speed).

Here $\rho$ is density and $c_s$ is the sound speed. Air gives $z_\text{air} = 1.2 \times 343 \approx 410$; water gives $z_\text{water} = 1000 \times 1481 \approx 1.5\times10^6$. A factor of 3600 apart.

Impedance measures how readily sound loads into a medium. When two impedances are close, the wave crosses smoothly. When they differ wildly, most of it bounces back. That jump at the water surface is why the shout reflects.

Put the rounded values back into $R=(z_2-z_1)/(z_1+z_2)$. With $z_1\approx 410$ and $z_2\approx 1.5\times 10^6$ the ratio is about 3600, $R$ sits on 0.999, and the square is the 99.9% above. $T=2z_2/(z_1+z_2)$ sits just under 2.

Reflection and Transmission — Decided by the Impedance Jump#

For a plane wave at normal incidence, the reflection coefficient $R$ and transmission coefficient $T$ depend only on impedance (paper eq. 37).

$z_1$ and $z_2$ are the incident-side and transmitted-side impedances. For air-to-water, $z_2 \gg z_1$, so $R \to 1$: almost total reflection. The pressure transmission coefficient $T$ can exceed 1. A larger pressure amplitude does not break energy conservation. The transmitted particle velocity drops, so the energy flux still balances.

The tube above already pushed that jump. A larger speed ratio and density ratio shrink the transmitted pulse and grow the reflected one. At 1:1 the reflection disappears. A different place where the same wave must not come back from an outlet is an NSCBC outflow.

The 13-Degree Wall — Snell's Law and the Critical Angle#

At oblique incidence the transmitted wave bends. The same Snell's law (the ratio of the sines of the incidence and transmission angles equals the ratio of sound speeds) that governs light holds for sound (paper eq. 41).

$\theta_i$ is the incidence angle and $\theta_t$ the transmission angle. When the transmitted medium is faster ($c_{s2} > c_{s1}$, like air to water), $\sin\theta_t$ exceeds 1 once the incidence angle passes a threshold. No transmitted wave can exist. That threshold is the critical angle.

$343/1481$ is about 0.232. The arcsine of that is 13.4 degrees. If an airborne source hits the surface more obliquely than this, the sound does not enter the water. Total reflection.

Move the incidence angle and the speed ratio, and the transmitted ray lies down on the surface or stands back up.

Push the incidence angle toward 13 degrees and the transmitted ray lies almost flat before it disappears. Only the reflected ray remains. Lower the speed ratio and the critical angle grows, so a steeper incidence still transmits.

Loading Sound onto a Weak-Compressibility Model#

An incompressible solver uses pressure only as a tool to enforce the divergence constraint. That pressure field cannot carry sound. The paper introduces a weak-compressibility (a time-derivative term added to pressure that grants a finite sound speed) equation.

$p$ is pressure, $\mathbf{u}$ is velocity, and $c_s$ is the local sound speed. The edited spot is the per-phase interpolation of $c_s$. Density and sound speed mix through the concentration $c$.

With that change, the same solver carries sound at $c_{s1}$ in medium 1 and $c_{s2}$ in medium 2. The authors discretize it with an entropy-stable discontinuous Galerkin spectral element method (DGSEM) and show spectral convergence.

The Modeling Error a Diffuse Interface Introduces#

Unlike VOF or level set, which cut the interface into a single cell, the Cahn–Hilliard model smears it over a width $\varepsilon$.

This smearing introduces a modeling error. The paper confirms numerically that the error decays as the square of the interface width — precisely $\mathcal{O}\!\left((\varepsilon/\lambda)^2\right)$ once normalized by wavelength. As $\varepsilon \to 0$, it converges to the sharp-interface analytic solution. The "interface width" slider on the 1D tube is this $\varepsilon$. Crank it up and the transmitted and reflected pulses smear slightly.

Python — Reproducing the Coefficients in a 1D Tube#

The linear acoustic equations are integrated with a staggered-grid FDTD scheme. A Gaussian pulse is fired, and the reflected and transmitted amplitudes after the interface are compared with the analytic values.

import numpy as np
 
def impedance(rho, c):
    """Acoustic impedance z = rho * c (density x sound speed)"""
    return rho * c
 
def analytic_coeffs(z1, z2):
    """Reflection / transmission for a normal-incidence plane wave (paper eq. 37)"""
    R = (z2 - z1) / (z1 + z2)
    T = 2.0 * z2 / (z1 + z2)
    return R, T
 
def tanh_interface(x, eps):
    """Concentration field smearing the interface over width eps (paper eq. 32)"""
    return 1.0 - 0.5 * (1.0 + np.tanh(2.0 * x / eps))
 
def run_acoustic_tube(rho1=1.0, c1=1.0, rho2=2.5, c2=1.8,
                      N=2000, L=4.0, eps=0.04, steps=2600):
    """Integrate the 1D linear acoustic system and measure pulse amplitudes."""
    dx = L / N
    x = -L / 2 + (np.arange(N) + 0.5) * dx
    conc = tanh_interface(x, eps)            # 1: medium 1, 0: medium 2
    rho = rho1 * conc + rho2 * (1.0 - conc)  # density interpolation
    cs  = c1   * conc + c2   * (1.0 - conc)  # sound speed interpolation
    K   = rho * cs**2                        # bulk modulus rho c^2
 
    p = np.exp(-((x + 0.55) / 0.12)**2)      # right-going Gaussian pulse
    u = p / (rho1 * c1)                      # medium-1 characteristic: u = p/(rho c)
 
    dt = 0.4 * dx / cs.max()                 # CFL-limited time step
    left, right = N // 4, 3 * N // 4         # left / right probes
    p_inc = p.max()
    refl_peak = tran_peak = 0.0
    for n in range(steps):
        u[:-1] += -dt / (0.5 * (rho[:-1] + rho[1:]) * dx) * (p[1:] - p[:-1])
        p[1:]  += -dt * K[1:] / dx * (u[1:] - u[:-1])
        if n > steps // 2:                   # after the pulse crosses
            refl_peak = max(refl_peak, np.abs(p[:left]).max())
            tran_peak = max(tran_peak, np.abs(p[right:]).max())
    return p_inc, refl_peak, tran_peak
 
if __name__ == "__main__":
    rho1, c1, rho2, c2 = 1.0, 1.0, 2.5, 1.8
    z1, z2 = impedance(rho1, c1), impedance(rho2, c2)
    R, T = analytic_coeffs(z1, z2)
    p_inc, refl, tran = run_acoustic_tube(rho1, c1, rho2, c2)
    print(f"z1={z1:.3f}  z2={z2:.3f}")
    print(f"analytic   |R|={abs(R):.3f}  T={T:.3f}")
    print(f"measured   |R|={refl / p_inc:.3f}  T={tran / p_inc:.3f}")

At $z_1=1.0$ and $z_2=4.5$, the measured values track the analytic $|R|=0.636$ and $T=1.636$. The defaults are $\rho_1=1$, $c_1=1$, $\rho_2=2.5$, $c_2=1.8$, so $z_2=\rho_2 c_2=4.5$. Then $R=(4.5-1)/(4.5+1)=3.5/5.5=0.636$ and $T=2\times 4.5/5.5=1.636$. The square of $0.636$ is near 0.405. $T=1.636$ does not push reflected energy past 1. Particle velocity drops with the impedance ratio and the flux still balances.

The grid is $N=2000$, $L=4$, so $dx=0.002$. $\varepsilon=0.04$ spans 20 cells, above the 10 points packed into the 3D air–water interface. With $c_2=1.8$, $dt=0.4\,dx/1.8$. Peaks are recorded after step 1300 of 2600.

What Felt Shaky While Reproducing It#

The first trap was separating the reflected wave. The left probe sees the sum of incident and reflected pressure. The paper runs a separate single-phase simulation to subtract the incident part. The code above measures only the residual on the left after the pulse has left, but it is off while the incident and reflected waves overlap.

Weak compressibility is not true compressibility. It cannot handle shocks or high Mach numbers. The authors restrict it to low-Mach, linear acoustics. The purpose differs from solving acoustics directly with a compressible solver.

The 3D air-water example used 67.5 million degrees of freedom, because the diffuse interface needs more than ten points packed into it. The cost explodes as the interface thins. Whether this replaces VOF depends on the problem scale.

A Paper to Read Next#

Shrinking $\varepsilon$ cuts the error as $(\varepsilon/\lambda)^2$. In 3D the demand of 10 points inside the interface pushed the count to 67.5 million degrees of freedom. The 1D default $\varepsilon/dx=20$ is twice that 10-point floor. Before cutting thinner, check that 0.636 and 1.636 already match. If they do not, the probe time is the problem, not anti-diffusion interface sharpening.

The same group's entropy-stable DGSEM iNS/CH paper (Manzanero et al., 2020) is the discretization this weak-compressibility term sits on. An outlet that swallows sound meets the NSCBC note above. Before forcing the interface thinner, check that normal incidence at an impedance ratio of 3600 puts $R$ on 0.999.