A common approach to account for current effects is to use an effective sound speed ceff, which is to add the current directly to the sound speed profile, which interprets current effects in terms of an effective sound speed profile (ESSP). This approach is often used in the literatures [25, 26].
The non-rigorous nature of introducing ceff is because the refractive index of a moving medium is not as simple as stationary medium. When introducing effective sound speed, the real moving medium is often replaced by a hypothetical station-ary medium with sound speed equal to effective sound speed. By this approach, the general ray equations in a stationary medium is applicable, which simplify the treatment of ray-tracing in a moving medium. And because the typical current ve-locity is much less than sound speed in ocean, the sound propagation is dominate by sound speed, this approach seems to be persuasive.
Generally, the effective sound speed ceff is defined as [20]
ceff = c + vR, (2.32)
where c is the real sound speed, vR= v· n is the component of the current velocity along the direction of propagation.
For some cases, e.g. the direction of sound propagation coincides with the direction of current. A further simplification can be made that
cs-eff = c + v, (2.33)
where the current magnitude is added directly to the sound speed, without projecting to n.
Figure 2.4: (a) Sound speed profile (b) Current profile (c) Simplified effective sound speed profile.
However, it is not physically correct to account for current effect using the effective sound speed approach. Figure 2.5(a) shows the ray velocity vector in a moving medium and Figure 2.5(b) shows how the effective sound speed approach accounts for the ray velocity. As depicted in Figure 2.5(a), the ray velocity in a moving medium is cn + v, the sum of two vectors. Whereas in the effective sound speed approach, the ray velocity is approximated by c + v · n. There are two differences between the effective sound speed approach and a correct ray-tracing for a moving medium. First, the wavefront normal vector n and the actual ray path cn + v are not collinear, which means that the effective sound speed cannot account for the shifting of wavefront by current. Second, the magnitude of ray velocity in a moving medium|cn + v| is different from the effective sound speed c + v · n.
It can be shown, via the Cauchy-Schwarz and triangular inequalities, that the relationship between average sound speeds calculated by correctly including the current effects in the ray tracings and those by two simplified approaches is
c + v· n
� �� �
ceff
≤ |cn + v| ≤ c + v� �� �
cs-eff
. (2.34)
Equation (2.34) shows that the effective sound speed approach will under-estimate the actual average speed along ray path, whereas the simplified effective sound speed approach will over-estimate it.
xp v
cn
cn + v
cn
Horizontal xp
v v · n
c
(a) (b)
Figure 2.5: Schematic plots of the ray velocity vector modeled by (a) moving-medium ray-tracing and (b) the effective sound speed approach.
Introducing effective sound speed has greatly simplified the treatment of ray-tracing in a moving medium. However, the effective sound speed approach is not applicable for all conditions. Ostashev [20] had discussed for the conditions that the rigorous introduction of ceff is possible. It is valid only for short-range propagation.
Section 4.4 shows the numerical results using effective sound speed approach and discusses the modeling errors in ray proeperties such as ray length, average speed along ray and travel time for a moving medium.
Chapter 3
Numerical Simulation
This chapter describes the details of numerical simulation, including the Cur-rent Ray-Tracing program, the environments used in the simulation, the validation of the Current Ray-Tracing program and perturbations of bathymetry.
3.1 Current Ray-Tracing Program
This section states the details of Current Ray-Tracing program. The computer program is written in double precision using matlab. The program is based on the ray equations for a moving medium, Equations (2.20)–(2.23) [Section 2.2]. Given the source position and launching angle θ0 as the initial conditions, the ray trajectory and travel time can be obtained by integrating Equations (2.20)–(2.23). A higher-order Runge-Kutta method with variable step size (ode45 solver) is used to perform the range integration. The derivatives of depth, wavefront normal and travel time with respect to range are calculated at each step.
3.1.1 Integration of the Differential Equations
The ode45 solver [10] in matlab is used to perform the integration of Equa-tions (2.20)–(2.23). The ode45 solver is based on an explicit Runge-Kutta(4,5) formula, the Dormand-Prince pair [4], which combines a fourth order method and a fifth order method. Given a differential equation written as:
y� = f (x, y), (3.1)
y(xn) = yn. (3.2)
To predict a new value yn+1 by a step size h, we first compute four auxiliary quan-tities k1, k2, k3, k4 and use them to compute yn+1
k1 = hf (xn, yn), (3.3)
k2 = hf (xn+ h
2, yn+k1
2), (3.4)
k3 = hf (xn+ h
2, yn+k2
2), (3.5)
k4 = hf (xn+ h, yn+ k3), (3.6) yn+1 = yn+ 1
6(k1+ 2k2+ 2k3+ k4). (3.7) The modified Runge-Kutta varies the step size h at each step to ensure the error is small enough. The error tolerance |ek| is controlled by two parameters: relative tolerance RelT ol and absolute tolerance AbsT ol, as follows:
|ek| ≤ max(RelT ol × |yk|, AbsT ol) (3.8) The relative and absolute tolerances are set as 1× 10−9 in the program for the consideration of accuracy and computational efficiency.
3.1.2 Sound Speed and Horizontal Current Velocity
The sound speed and horizontal current velocity of the waveguide vary in depth and range. These environmental parameters are given in gridded data. The values of sound speed and horizontal current velocity within a grid cell are obtained using cubic spline interpolation. The differentiation terms �∂c
∂r, ∂z∂c, ∂v∂r, & ∂v∂z�
in Equa-tions (2.21) and (2.22) are evaluated using central difference.
3.1.3 Boundaries and Reflection
The top of the medium is assumed to be pressure-release boundary, and the bottom is assumed to be hard boundary. The sea surface is assumed to be flat. For an irregular bottom, the bathymetry is given in a regular grid. The depth of an irregular bottom within a grid cell is calculated by linear interpolation.
At each step the program detects if a boundary reflection occurs in the next range step. If an impact has occurred, the initial ray position, travel time and θ for the next step are modified accordingly. The initial ray position and travel time are the values when the previous ray section hits the boundary. As for the initial θ, it is calculated using the reflection law to change the sign of θ after reflection. For a
flat bottom, the sin θ and cos θ after reflection are modified as follows
sin θr = − sin θi (3.9)
cos θr = cos θi (3.10)
where θr is the new initial θ for the next step, and θi is θ before reflection. For the case of irregular bottom, Figure 3.1, according to the law of reflection, θi−α = θr+α, the sin θ and cos θ after reflection are modified as follows
sin θr = − sin(θi− 2α) (3.11)
cos θr = cos(θi− 2α) (3.12)
where α is the angle of the bottom with respect to the horizontal. The slope of the interface is interpolated from the neighboring grid points of the impact point.
r i
Boundary normal
Horizontal Incident ray
5HÁHFWHGUD\
Bottom
surface
Figure 3.1: Schematic of ray reflecting from a piecewise linear boundary