Propagation in the Ocean Waveguide

Now we are going to include the effects of the seafloor. We will again assume that the sound speed is constant throughout the water column and that the seafloor and sea surface are both perfectly smooth. We will use the fluid approximation for the seafloor and assume that it has the properties of one of the three types of sediments we considered earlier. In order to focus just on the impacts of the sea surface and seafloor we will also neglect absorption in the water column. (It would be nice at a future date to provide the option to include the absorption as a function of frequency so that the user can see that effect as well.) In this example, each path between the source and receiver is referred to as a ray and script allows us to increase or decrease the number of rays included in the calculation. If the number of rays is set to 2, then the script will only include the direct and surface paths making it identical to the seafloor-only case considered above.

Show code cell source

Hide code cell source

# Import packages
import numpy as np
import matplotlib.pyplot as plt
from math import ceil, log, pi, atan, sqrt
from scipy.signal.windows import tukey
import ipywidgets as widgets

# Set inline
%matplotlib inline

Show code cell source

Hide code cell source

def CW_pulse_propagation(p, t, tv, frequency, tau, win):
    """
    Function to compute continous wave pulse using Tukey window.
    """
    # Compute the full pulse with the phase factor:
    CW_pulse = p * np.exp(-1j * (2*np.pi*frequency*t - np.pi/2))
    
    # Identify indices where the pulse is active: (t - tv) in [0, tau]
    ind = (t - tv >= 0) & (t - tv <= tau)
    
    # Zero out values outside the active pulse window
    CW_pulse[(t - tv) < 0] = 0
    CW_pulse[(t - tv) > tau] = 0
    
    if np.any(ind):
        # Apply the Tukey window
        window = tukey(np.sum(ind), win)
        CW_pulse[ind] *= window
    
    return CW_pulse

Show code cell source

Hide code cell source

def RCoeff(frequency, theta_i_degrees, cw, cp1, rho_ratio, rhow, alphap1):
    """
    Function to compute reflection coefficient.
    """
    # Convert attenuation to loss parameter
    deltap1 = alphap1 * cp1 * np.log(10) / (40 * np.pi * frequency)

    # Convert from degrees to radians
    theta_i = np.deg2rad(theta_i_degrees)

    # Complex sound speed
    cpp = cp1 / (1 + 1j * deltap1)

    # Ratio of sound speeds
    ap = cpp / cw

    # Impedance for water
    zw = rhow * cw / (np.abs(np.sin(theta_i)) + 1e-6)

    # Impedance of sediment
    zp = rho_ratio * cpp / np.sqrt(1 - (ap * np.cos(theta_i))**2)

    # Impedance ratio
    zwp = zp / zw

    # Reflection coefficient
    R1 = (zwp - 1) / (zwp + 1)

    return R1

Show code cell source

Hide code cell source

def simulate_multiple_ray_propagation(
    sediment_type,
    seafloor_depth,
    source_depth,
    receiver_depth,
    receiver_range,
    number_rays,
    time_range,
):
    """
    Simulates multiple ray propagation.
    """
    # Get start and end display times
    display_time_start, display_time_end = time_range

    # Frequency (Hz)
    f = 1000

    # Source level (dB)
    source_level = 0

    # Where the source starts
    source_range = 0

    # CW pulse cycles
    pulse_width = 8

    # Plot start and end times
    display_time_start = display_time_start
    display_time_end = display_time_end
    
    # Set propagation parameters based on sediment type
    if sediment_type == "Fine Sand":
        cw = 1500
        cp = 1.15 * cw
        rhop = 2.051
        rhow = 1
        deltap = 0.01
    elif sediment_type == "Silty Clay":
        cw = 1500
        cp = 0.982 * cw
        rhop = 1.546
        rhow = 1
        deltap = 0.001
    elif sediment_type == "Clayey Silt":
        cw = 1500
        cp = 1.01 * cw
        rhop = 1.597
        rhow = 1
        deltap = 0.005
    
    # dB/m conversion parameter
    alphap = 40 * pi * f * deltap / (cp * log(10))

    # Wavenumber in water
    k = 2 * pi * f / cw

    # Calculate rays (in multiples of 4)
    number_raysets = int(ceil(number_rays / 4))

    # Travel times for rays
    t_v_all = []

    # Pressures for rays
    p_v_all = []

    # Ray path coordinates
    rays = []

    # Iterate through ray sets
    for j_r in range(1, number_raysets + 1):
        # Calculate depth differences for virtual sources
        z_v = np.zeros(4)
        z_v[0] = 2 * seafloor_depth * (j_r - 1) + source_depth - receiver_depth
        z_v[2] = 2 * seafloor_depth * j_r - source_depth - receiver_depth
        z_v[1] = -(2 * seafloor_depth * (j_r - 1) + source_depth + receiver_depth)
        z_v[3] = -(2 * seafloor_depth * j_r - source_depth + receiver_depth)
        t_v = np.zeros(4)
        p_v = np.zeros(4, dtype=complex)
        r_v = np.zeros(4)
        r_angle = np.zeros(4)
        rc1 = np.zeros(4, dtype=complex)
        
        for j_v in range(4):
            # Grazing angle
            r_angle[j_v] = atan(abs(z_v[j_v]) / (receiver_range - source_range))

            # Reflection coefficient
            rc1[j_v] = RCoeff(f, np.degrees(r_angle[j_v]), cw, cp, rhop / rhow, rhow, alphap)

            # Ray travel distance
            r_v[j_v] = sqrt((receiver_range - source_range)**2 + z_v[j_v]**2)

            # Travel time
            t_v[j_v] = r_v[j_v] / cw
        
        # Calculate pressure for each ray
        p_v[0] = ((-rc1[0])**(j_r - 1)) * np.exp(1j * k * r_v[0]) / r_v[0]
        p_v[1] = -((-rc1[1])**(j_r - 1)) * np.exp(1j * k * r_v[1]) / r_v[1]
        p_v[2] = ((-rc1[2])**(j_r - 1)) * rc1[2] * np.exp(1j * k * r_v[2]) / r_v[2]
        p_v[3] = -((-rc1[3])**(j_r - 1)) * rc1[3] * np.exp(1j * k * r_v[3]) / r_v[3]
        
        # Determine intersections with the surface/seafloor

        # Slope of the ray line
        m_v = -z_v / (receiver_range - source_range)

        # Intercept
        b_v = receiver_depth - m_v * receiver_range
        
        # Calculate intersection levels
        l1 = np.arange(2 * (j_r - 1), 0, -1) * seafloor_depth
        l3 = np.arange(2 * j_r - 1, 0, -1) * seafloor_depth
        if j_r - 1 > 0:
            l2 = -np.append(np.arange(2 * (j_r - 1), 0, -1) * seafloor_depth, 0)
        else:
            l2 = np.array([0])
        l4 = -np.append(np.arange(2 * j_r - 1, 0, -1) * seafloor_depth, 0)
        
        # Compute intersections if the slope is nonzero
        if m_v[0] != 0 and l1.size > 0:
            r_i1 = (l1 - b_v[0]) / m_v[0]
            z_i1 = np.tile(np.array([0, seafloor_depth]), j_r - 1) if j_r - 1 > 0 else np.array([])
        else:
            r_i1, z_i1 = np.array([]), np.array([])

        if m_v[1] != 0 and l2.size > 0:
            r_i2 = (l2 - b_v[1]) / m_v[1]
            z_i2 = np.concatenate(([0], np.tile(np.array([seafloor_depth, 0]), j_r - 1)))
        else:
            r_i2, z_i2 = np.array([]), np.array([])

        if m_v[2] != 0 and l3.size > 0:
            r_i3 = (l3 - b_v[2]) / m_v[2]
            z_i3 = np.concatenate(([seafloor_depth], np.tile(np.array([0, seafloor_depth]), j_r - 1)))
        else:
            r_i3, z_i3 = np.array([]), np.array([])

        if m_v[3] != 0 and l4.size > 0:
            r_i4 = (l4 - b_v[3]) / m_v[3]
            z_i4 = np.tile(np.array([seafloor_depth, 0]), j_r)
        else:
            r_i4, z_i4 = np.array([]), np.array([])

        # Save ray paths
        rays.append({"ray_r": np.concatenate(([source_range], r_i1, [receiver_range])),
                     "ray_z": np.concatenate(([source_depth], z_i1, [receiver_depth]))})
        rays.append({"ray_r": np.concatenate(([source_range], r_i2, [receiver_range])),
                     "ray_z": np.concatenate(([source_depth], z_i2, [receiver_depth]))})
        rays.append({"ray_r": np.concatenate(([source_range], r_i3, [receiver_range])),
                     "ray_z": np.concatenate(([source_depth], z_i3, [receiver_depth]))})
        rays.append({"ray_r": np.concatenate(([source_range], r_i4, [receiver_range])),
                     "ray_z": np.concatenate(([source_depth], z_i4, [receiver_depth]))})
        
        # Save travel times and pressures
        t_v_all.extend(t_v.tolist())
        p_v_all.extend(p_v.tolist())
    
    # Time series calculation

    # Time step (T/10)
    dt = (1 / f) / 10
    t1 = np.arange(0, max(t_v_all) * 1.5, dt)
    tau = pulse_width / f      

    # Create CW pulse with a 20% Tukey window taper
    output_pulse = np.zeros_like(t1, dtype=complex)
    for j in range(min(number_rays, len(p_v_all))):
        output_pulse += CW_pulse_propagation(p_v_all[j], t1, t_v_all[j], f, tau, 0.2)
    
    # Create figure containing three subplots
    _, axes = plt.subplots(1, 3, figsize=(18, 5))

    # Create ray paths plot
    for j in range(min(number_rays, len(rays))):
        axes[0].plot(rays[j]["ray_r"], rays[j]["ray_z"], "k")
    axes[0].invert_yaxis()
    axes[0].set_xlim(source_range, receiver_range)
    axes[0].set_ylim(0, seafloor_depth)
    axes[0].set_xlabel("Range (m)")
    axes[0].set_ylabel("Depth (m)")
    axes[0].set_title("Ray Paths")
    axes[0].grid(True)

    # Create pressure time series plot
    axes[1].plot(t1, np.real(output_pulse) * np.sqrt(2) * 10**(source_level / 20), linewidth=0.3)
    axes[1].set_ylim(-0.01, 0.01)
    axes[1].set_xlim(display_time_start, display_time_end)
    axes[1].set_xlabel("Time (s)")
    axes[1].set_ylabel("Pressure (µPa)")
    axes[1].set_title("Pressure Time Series")
    axes[1].grid(True)

    # Create intensity time series plot
    intensity = 20 * np.log10(np.abs(output_pulse) + 1e-6) + source_level
    axes[2].plot(t1, intensity, linewidth=2)
    axes[2].set_ylim(-65, -30)
    axes[2].set_xlim(display_time_start, display_time_end)
    axes[2].set_xlabel("Time (s)")
    axes[2].set_ylabel("Intensity (dB)")
    axes[2].set_title("Intensity Time Series")
    axes[2].grid(True)

    # Show plots
    plt.tight_layout()
    plt.show()


# Create interactive widgets
layout = widgets.Layout(width="500px", height="30px")
_ = widgets.interact(
    simulate_multiple_ray_propagation,
    sediment_type=widgets.Dropdown(
        options=["Fine Sand", "Silty Clay", "Clayey Silt"],
        value="Clayey Silt",
        description="Sediment Type",
        style={"description_width": "initial"}, 
        layout=layout,
    ),
    seafloor_depth=widgets.IntSlider(
        min=30,
        max=100,
        step=1,
        value=50,
        description="Seafloor Depth (m)",
        style={"description_width": "initial"}, 
        layout=layout,
    ),
    source_depth=widgets.IntSlider(
        min=0,
        max=50,
        step=1,
        value=25,
        description="Source Depth (m)",
        style={"description_width": "initial"}, 
        layout=layout,
    ),
    receiver_depth=widgets.IntSlider(
        min=0,
        max=100,
        step=1,
        value=20,
        description="Receiver Depth (m)",
        style={"description_width": "initial"},
        layout=layout,
    ),
    receiver_range=widgets.IntSlider(
        min=1,
        max=500,
        step=10,
        value=100,
        description="Receiver Range (m)",
        style={"description_width": "initial"},
        layout=layout,
    ),
    number_rays=widgets.IntSlider(
        min=4,
        max=100,
        step=4,
        value=10,
        description="Number of Rays",
        style={"description_width": "initial"},
        layout=layout,
    ),
    time_range=widgets.FloatRangeSlider(
        value=[0, 0.25],
        min=0,
        max=5,
        step=0.01,
        description="Time Range (s):",
        style={"description_width": "initial"},
        layout=layout,
    )
)

As more rays are added to the calculation, the angles of the new rays get steeper and steeper. This will cause the seafloor reflection coefficient to decrease, increase the number of times the ray will reflect from the seafloor, and will meant that each path will get longer and longer. All of these will lead to a decrease in the pressure level for the sound traveling along that path. As a consequence, the contributions from additional rays will become less and less important to the received signal. To see this, we will use the ray paths to calculate the time series for the source transmitting a single frequency pulse, similar to what was done in the seafloor-only case.

As before, if the pulse is short, high frequency, and near the source, the individual arrivals coming in along the different paths can be clearly discerned. As the frequency is decreased and the pulse is lengthened, the individual paths interfere. If enough rays are added, the contributes smooth out to a steady level.

A more realistic ocean environment could include a depth-dependent sound speed in the water column, range-dependence to the bathymetry of the seafloor, changes in the seafloor properties as a function of range and depth, etc. With a more complicated environment comes the need for more involved computational methods to model the acoustic field. Future notebooks will explore these methods and provide examples.