Ray dynamics 69 multiple sectioning of the trajectories yielded figures similar to Poincar´e sections of pe-riodically perturbed trajectories (Figure 4.6). The procedure is detailed at the end of the previous section. It is pointed out here that given a certain range of simulation, the number of points to plot is increasing with increasing δ. The accuracy, however, whether these points in phase space are close to the sectioning surface is decreasing. In contrast, with the approach presented above, the points which constitute the magnified curve in Figure 4.13 lie on the surface of intersection at the axis depth to within the accuracy of a root-finding procedure, and the range of simulation solely determines the number of points to plot. Furthermore, Figure 4.13 supports the idea of nonintersect-ing regular and chaotic regions in the phase space of the ray equations, i.e. a dense set of stable solutions. Note the cascade of sections of regular orbits.
It is pointed out that increasing N implies higher dimensionality of the problem.
Brown (1998) associated this with the higher [2(N+ 1)] dimensional phase space of the equivalent autonomous system. Here it is manifested by the following. We consider regular trajectories in the representation of the third reduction technique in case of N
= 0, 1, 2. With no perturbation (N = 0), sectioning the trajectories results in single points; with harmonic perturbation (N = 1) the sections of trajectories are curve segments. In case of the simplest quasiperiodic perturbation (N = 2) the increasing dimensionality is indicated by the intersecting sections of regular trajectories (Figure 4.13).
To analyse the divergence of nearby trajectories two different versions of the varia-tional equations can be derived respectively for the deformation gradient X and the perturbationδx. Equation (4.97) is differentiated with respect tox0 which implies that
∂
∂x0 µ∂x
∂r
¶
= ∂f(x, r)
∂x · ∂x
∂x0. (4.99)
Differentiation with respect to variablesr and x0 is interchangeable due to their inde-pendence. Utilising this results in the variational equation written for the deformation gradient:
X0 = ∂f(x, r)
∂x ·X. (4.100)
The ‘initial’ condition for the variational equation yields in form of the identity matrix by differentiating the initial condition of the original system. Differentiation of the left hand side yields the new definiendum: X(x0,0) = ∂x(x0,0)/∂x0, and differentiation of the right hand side yields the new definien: ∂x0/∂x0 = I, hence, X(x0,0) = I.
A small perturbation of the dependent variables is introduced in tune with complete differentials as:
δx = ∂x
∂x0δx0 =X·δx0. (4.101)
It is worth mentioning that the above expression constitutes a linear mapping between initial and actual values of the perturbation. Differentiating with respect to r, and utilising (4.100), one obtains a system of differential equations for the perturbation too:
δx0 = ∂f(x, r)
∂x ·δx. (4.102)
Equation (4.100) is to define a measure of stability, the Lyapunov exponents; (4.102) will be applied to a simplified calculation. They are both linear systems of differential equations, with an identical nonconstant coefficient matrix. First, local stability is analysed by the coefficient matrix.
At any range, a system with some constant coefficient matrix A can be associated with (4.102). When the coefficient matrixA is obtained by evaluating ∂f(x, r)/∂x at that particular range, the two systems are qualitatively equivalent at this point locally.
Hence, the local stability of the original system can be analysed by the associated one.
For such a system the solution is known in a closed exponential form: Σni=1ciuiexp(λir), whereλi and ui are respectively eigenvalues and eigenvectors of the coefficient matrix A, and ci are constants determined by initial condition. The eigenvalues determine
Ray dynamics 71 stability, such that the solution is unstable when there is at least one eigenvalue with positive real part. (Weisstein n.d. First-Order Ordinary Differential Equation, online content)(Arfken 1985, p. 440-445)
The coefficient matrix in case of the above introduced ray system of Hamiltonian (4.7) is found to be:
"
0 c0
−c10∂∂z2V2 0
#
with eigenvalues λ1,2 = ± q
−∂∂z2V2. The eigenvalues come in either a plus-minus real pair, or as complex conjugates. For the latter, a small neighbourhood of phase space volume is rotating around the reference trajectory in the phase space of ray trajectories;
for the former, the phase space volume is stretching in one direction and shrinking in another one. By the exponential form of the solution, the real positive eigenvalue is regarded as the Exponential Increment of Local Instability (EILI); with complex conjugate eigenvalues this is zero. Thus,
EILI : Re r
−∂2V
∂z2 . (4.103)
Yan (1993) arrived at an equivalent form of the criterion for local instability.
Long range instability, on the other hand, is analysed by the deformation gradient.
The exponential form Σni=1ciuiexp(λir) which locally approximates the perturbation δx is derived from the general solution which is in form of a matrix exponential: eAr. According to the eigen decomposition theorem, eigenvalues of the matrix exponential are exponentials of the eigenvalues of matrixA (Weisstein n.d. Eigen Decomposition, online content)(Press et al. 1992, p. 51-63):
λX,i=eλir. (4.104)
Since ∂xf(x, r) generally varies with range, eigenvalues of the deformation gradient do not have such a simple exponential form. However, the above representation is suitable to motivate the following formula for the average exponential growth of its singular values, called Lyapunov exponents:
νi = lim
r→0
ln(si)
r . (4.105)
The singular valuessi are the square roots of the real eigenvalues of the matrixX·XT, which are always nonnegative. They give the measure of the principal axes of the error ellipsoid associated with the stretching and shrinking of phase space volume. Note
0 100 200 300 400 0
0.1 0.2
EILI[1/km]
0 100 200 300 400 0 100 200 300 400
0 2000 4000 6000 0
0.05 0.1 0.15 0.2
νL[1/10km]
2000 4000 6000 range [km]
2000 4000 6000
Figure 4.14: Exponential increment of local instability (above) and asymptotic eval-uation of Maximal Lyapunov Exponent (bellow) of three rays: on the sides, the two rays from Fig. 4.1, in between, a ray in case of the double duct profile with initial conditions (z0, ϕ0) = (0.325 km,−7◦) already cited in Fig. 4.2, as well as the other two rays
also the following property of the singular values. Recall that the the Jacobian of the one-degree-of-freedom Hamiltonian flow is constantly unity due to its initial value:
J0(x0,0) = det(I) = 1, which follows from (4.23) and (4.24). Thus, det(X ·XT) = det(X)det(XT) = det2(X) = 1 too. That is, s21s22 = 1, and also s1s2 = 1. Therefore, ln(s1) = −ln(s2). The last equality suggests that the Lyapunov exponents come in a plus-minus pair, i.e. ν1 =−ν2. That is, the larger is either zero or positive, denoting stable or chaotic motion.
In practice, it is not feasible to calculate the Lyapunov exponents by the above formula (4.105), but different numerical algorithms may be employed. The algorithm applied to underwater acoustics by Smithet al. (1992b) will now be discussed.
Equations (4.102) along with the original equations (4.97) are integrated to find the perturbation, by which the distance between reference and perturbed trajectories is approximated by the following measure: d=|δz|+|δp|. By the latter definition it is clear that the variational equation must be nondimensional. These equations for the ray system are taken to be:
Ray dynamics 73
Figure 4.15: Maps of the MLE for the Munk (left) and double duct profiles (right), with resolution given by a grid of 200×200 initial conditions in a regime: 0 < z0 <4 km and−15◦ < ϕ0 <15◦
δz0 = δp, (4.106)
δp0 = −R2∂z,zV(r, z)δz, (4.107) where R, the internal wave length was chosen as normalisation constant for creating nondimensional variables fromrand z. Initial conditions are chosen arbitrarily, but so thatd(0) = 1. The larger or Maximal Lyapunov Exponent (MLE) is then approximated by taking
νL = ln(d)
r (4.108)
at some sufficiently large r. The calculation might need to be reinitialized to avoid overflow problems. Note that the Gram-Schmidt technique (Alligood et al. 1997) treats the MLE separately from the others, and only the measure d is approximated above, which isp
δz2+δp2 originally. Nondimensionalization of depth and range was achieved dividing by the internal wave length; andp is already nondimensional.
Figure 4.14 illustrates the application of the short and long term concepts of stability via three examples. The left upper and lower diagrams correspond with a stable ray.
There is no local instability, and the MLE estimate apparently converges to zero, but very slowly. The case presented in the middle pair is an interesting one with periodic recurrent local instabilities, yet it is of a stable ray having zero MLE. The third pair of diagrams corresponds with a chaotic ray; aperiodic recurrent local instabilities go with a nonzero MLE. Furthermore, this pair illustrates the difficulties with the finite range estimation of the MLE as the asymptotic value is not clear. It is concluded thus that recurrence of local instabilities is not a sufficient condition for long term instability.
Note the regular pattern of recurrence in the stable case, as opposed to the apparently irregular pattern of the chaotic one.
Based on these calculations we can develop a method to explore the space of initial conditions with respect to stability by numerical means. This is done by estimating the MLE associated with a grid of initial conditions. Figure 4.15 presents results for models using the two sound speed profiles introduced earlier. The equations were integrated over a range of 6000 km again. Reflection of rays from a flat sea surface and sea bottom was taken into account by reversing the ray angle and the sign ofδp; the ranges of the points of reflection were determined exactly employing a trick devised by H´enon (1982). The integration was carried out using fourth order Runge-Kutta integrator with a fixed step size of 0.1 km.
Chapter 5
A more descriptive model
5.1 Internal wave induced sound speed perturba-tion
In the previous chapter the ray equations included a very simple driving term to model horizontal variability, harmonically varying with range and exponentially decaying with depth. It was interpreted as an approximation of a single mode of internal wave motion.
To depart from this highly idealised picture in the present chapter a physically realistic model of internal waves is introduced, and, as presented in (B´odai et al. 2009), an analysis of ray stability in single and double channel wave guides is pursued using this model.
Horizontal variability of the sound speed field can partly be accounted for the undulation of iso-density surfaces, a result of vertical water parcel displacement due to internal wave motion. We recall the fundamental law of ocean structure (2.7) which establishes this connection such as:
δc≈cµ
gN2ζ. (5.1)
In the above, ζ will be later defined as the vertical displacement: a perturbation quantity; and parameters appearing in (2.7) are lumped together in the nondimensional parameter µ (taken to be 24.5 in computations). Hence, in the following it will be sufficient to describe the kinematics of internal wave motion.
We start with the equations of motion for inviscid incompressible rotating flow (Flatt´e 1979, p. 47.)(Eriksen 1976, p. 53.):
ρDtv+fˆk×v=−p∇T +ρgˆk, (5.2)
Dtρ+ρ(∇ ·v) = 0, (5.3)
∇ ·v= 0, (5.4)
however variation in density is allowed for. Equation (5.2) is Euler’s equation in a rotating frame, wheref = 2Ω sin(lattitude) is the Coriolis parameter with Ω = 2π/day, the angular speed of the Earth, and ˆk is the unit vector directed upward from the surface of Earth.1 Our computations were carried out for a latitude of 30◦. Equation (5.3) expresses the law of continuity, and (5.4) is required for incompressibility. The symbol Dt( ) = ∂t( ) + v· ( )∇T defines the material time derivative. The depth coordinate z is positive downward, and let its 0 level coincide with the water surface.
A decomposition of the density and pressure into static and time-varying parts is proposed next, such that
p=p0+p0 and ρ=ρ0+ρ0, (5.5) where time-varying parts are denoted by primes, and the static parts are defined by the equations of hydrostatics:
p0∇T =ρ0gˆk. (5.6)
Utilising this and neglecting all second-order small terms, the original system of equa-tions reduces to the following linear system:
ρ0∂tu−f v =−∂xp0, (5.7)
ρ0∂tv+f u=−∂yp0, (5.8)
ρ0∂tw=−∂zp0+ρ0g, (5.9)
∂tρ0+w∂zρ0 = 0, (5.10)
∂xu+∂yv+∂zw= 0. (5.11)
The solution of the linear equations is sought in form of a plane wave, whereby all the variables have equal phase, and some depth variation, that is:
1In the vertical direction gravitational and pressure forces are dominating, and the vertical com-ponent of centrifugal forces is included in the gravitation; in horizontal directions the Coriolis forces generally dominate over centrifugal forces (Wells 1986, p. 194.)
A more descriptive model 77
w=W(z)ei(kxx+kyy−ωt), etc. and k2 =k2x+ky2. (5.12) Eliminating all terms butW in (5.7)-(5.11), and identifying g∂zρ0/ρ0 as the buoyancy frequency squaredn2(z), we have:
∂zzW +
µn2−ω2 ω2−f2
¶
k2W =−n2
g ∂zW. (5.13)
It can be shown (Flatt´e 1979, p. 48.) that the term on the right is always negligible.
According to measurements the longest vertical internal wave lengths are of 1000 m magnitude, which meanskV= 1 km−1 minimal vertical wave number. In this extremal case the ratio of (minus) the term on the right and one of the terms on the left is estimated as:
n2 g
∂zW
∂zzW ≈ n2
gkV <4×10−4. (5.14)
With this we may write the equation for the approximate modes of the vertical particle velocity as follows:
∂zzW +
µn2−ω2 ω2−f2
¶
k2W = 0, (5.15)
which, by integration, is also the equation for the modes of the vertical particle dis-placementζ.
To solve the modal equations, appropriate boundary conditions must be set. At the bottom with zero velocity this is W(H) = 0. It can be argued (Flatt´e 1979, p.
49.) that a similar condition is approximately true at the surface, such thatW(0) = 0.
With these the modal equation poses a classical Sturm-Liouville problem when only discrete values of ω and corresponding k lead to a solution, and so an infinite series of modes W(j, k, z) with mode numbers j = 1,2, . . . satisfy the equations. The modes are orthogonal so that
Z H
0
(n(z)2−f2)W(j, k, z)W(i, k, z)dz = 0 when i6=j. (5.16) As the modes are undefined up to a constant multiplier, an orthonormaility condition can be set as (Brown & Viechnicki 1998):
Z H
0
(n(z)2−f2)W(j, k, z)W(i, k, z)dz =δji. (5.17) The modal equation (5.15) has an identical form to Bessel’s equations when the buoy-ancy frequencyn(z) assumes an exponential form (2.11), and the solution can be found
in terms of Bessel functions of the first and second kind (Flatt´e 1979, p. 53.). With arbitrary stratification the modal equation can be solved numerically; or, for compu-tational efficiency, an approximate analytical solution can be employed. The one to be used in our computations, given by Colosi & Brown (1998), takes the following form:
Vj(ξ) = (1/n0)(2/B)1/2sin(jπξ) and kBn0/p
ω2−f2 =jπ, (5.18) where
V(z) = (n(z)/n0)1/2W(z), (5.19)
ξ(z) = 1 n0B
Z H
z
n(z0)dz0, and (5.20)
n0B = Z H
0
n(z)dz. (5.21)
The modes are normalised so that the orthonormaility condition (5.17) is approximately satisfied with n2 −f2 ≈ n2. Assuming exponential stratification, a condition of the approximation is found in the form: (2jπ)2 À(n0/n)2.
With the modes that form a complete set, a statistical description of the vertical displacement can be given:
hζ2i= 1 2
Z Z X
j
G(j, kx, ky)W2(j, k, z)dkxdky. (5.22) In the aboveG(j, kx, ky) represents the internal wave spectrum, first determined empir-ically by Garrett & Munk (1972,1975). The pair of angle bracketsh·idenotes average.
The errors that the approximate modes described above would introduce for small mode numbers can be compensated for by choosing appropriate spectral weights when matching experimental data. The internal wave spectrum consistent with the approx-imations of the vertical modes are given by Colosi & Brown (1998) as:
G(j, kx, ky) = 2B3n20E π2M
1 j2+j∗2
kjp
kx2+ky2
(k2x+ky2+k2j)2, (5.23) where kj = (f /n0)(πj/B), M = P∞
j=1(j2+j∗2)−1 ≈ 12j∗−2(πj∗ −1) is a normalisation constant, andj∗ = 3 and E = 6.3×10−5 are empirical dimensionless constants. With this, the vertical displacement is obtained explicitly in the following form:
ζ(x, y, z, t) = Re
"Z Z X
j
g(j, kx, ky)W(j, k, z)ei(kxx+kyy−ω(j,k)t)dkxdky
#
, (5.24)
A more descriptive model 79 whereg(j, kx, ky) is a complex Gaussian random variable with zero mean and variance:
hg(j, kx, ky)g∗(j0, kx0, ky0)i=G(j, kx, ky)δ(kx−k0x)δ(ky−ky0).
For our purposes, the spectrum is required in a vertical slice of the ocean and frozen in time. This can be obtained by integrating with respect to, say,ky, and substituting at, say, zero time. The arbitrary choice between kx and ky is partially justified by the horizontal isotropy already built into G(j, kx, ky) = G(j, k), which is valid if the rotation of the Earth is neglected. Taking a frozen slice is justified by pointing to the greatly different time scales of the internal wave and acoustic phenomena. This spectrum was originally derived by Brill & Dozier (1985):
G(j, kx) = Z ∞
−∞
G(j, kx, ky)dky = I(j, kx) j2+j∗2
2B3n20E
π2M , (5.25)
where
I(j, kr) = kj
kr2+kj2 +1 2
k2r
(k2r+kj2)3/2 ln
q
k2r+kj2+kj q
kr2+kj2−kj
. (5.26)
Finally we mention that Colosi & Brown (1998) pointed out that the formulation of the vertical displacement given by (5.24) suggests that its field can be numerically generated by a multidimensional (inverse) Fast Fourier Transformation. The formula we consider for the frozen vertical slice problem is assembled as:
ζ = Re
"
2B π
µE M
¶1/2³n0 n
´1/2
(∆kr)1/2X
j
sin(jπξ(z)) (j2+j∗2)1/2
X
kr
I1/2(j, kr)eiφ(j,kr)eikrr
# , (5.27) where ∆kr is the step size of the discretised formula of integration over the variable kr; φ(j, kr) are randomly chosen phase angles. The depth of the ocean is taken to be H = 4 km. To perform the numerical Fourier transformation, the variable of range and the horizontal wave number are discretised as rk = kh and kr,k = 2πk/(Nh), k = 0. . . N −1, where h = 0.1 km is the step size of the fourth order Runge-Kutta integrator applied to the ray system, and N is some power of 2. The upper limit of range yielding asNhis desired to exceed 5000 km, the range of ray tracing performed for our results. To obtain the minimal N now it has to be 6553.6 km, with which N = 216. (The MathWorks n.d. fft, Matlab Help article online) For our computations we used a maximum vertical mode number of 50. The minimal value ofkr here is set to be 2π/100 km. The (inverse) numerical Fourier transform was actually carried out by Matlab routineifft. With all the above indicated values the introduction of anad hoc
0
0.5
1
1.5
2
depth[km]
1.49 1.51 0.5
1
1.5
2
sound speed [km/s]
depth[km]
0 50 100 150 200 250
range [km]
−5 0 5
ray angle [◦]
Figure 5.1: Examples of ray paths (middle panel) and ray trajectories (left panel) with both Munk (above) and double duct profiles (below), shown on the left. The magenta and blue curves are respectively associated with the unperturbed and perturbed sys-tems. Initial conditions were: (z0, ϕ0) = (1 km,4◦) – Munk, (z0, ϕ0) = (1.3447 km,4◦) – double duct. Perturbation parameters as specified in the text – unchanged for all computations
division by 5 was necessary at (5.27) to achieve approximately the same amplitudes of the wave forms of the perturbation field as those presented by Colosi & Brown (1998).