Acoustic and geoacoustic inverse problems in randomly perturbed shallow-water environments
Laure Dumaz Josselin Garnier Guilhem Lepoultier






Ceremade, Université Paris-Dauphine,

75016 Paris - France,
and

Sivienn, 

29570 Roscanvel - France


Centre de Mathématiques Appliquées, Ecole Polytechnique, 91128 Palaiseau Cedex - France, 
and

Sivienn, 29570 Roscanvel - France



Sivienn, 29570 Roscanvel - France





The main goal of this paper is to estimate the regional acoustic and geoacoustic shallow-water environment 
from data collected by a vertical hydrophone array and transmitted by distant time-harmonic point sources.
We aim at estimating the statistical properties of the random fluctuations of the index of refraction
in the water column and the characteristics of the sea bottom.
We first explain from first principles how acoustic wave propagation can be expressed as Markovian dynamics for 
the complex mode amplitudes of the sound pressure, which makes it possible to express the cross moments of the sound pressure
in terms of the parameters to be estimated.
We then show how the estimation problem can be formulated as a nonlinear inverse problem 
using this formulation, that can be solved by minimization of a misfit function. 
We apply this method to experimental data collected by the ALMA system (Acoustic Laboratory for 
Marine Applications).





43.30.Bp, 
43.30.Re, 	
43.20.Mv, 
43.20.Bi. 




Acoustic and geoacoustic inverse problems in randomly perturbed shallow-water environments
    Guilhem Lepoultier
    December 30, 2023
==========================================================================================





§ INTRODUCTION



In this paper we consider acoustic wave propagation in a randomly perturbed
shallow-water waveguide with an absorbing sea bottom.
The random perturbations of the index of refraction are due to internal waves,
which induce  temperature  and  salinity  fluctuations, 
and the sea bottom  is made of sediments, which induce dissipation.
In the regime where random perturbations and dissipation are small and propagation distance is large,
it is possible to get an effective description of the acoustic wave propagation
in terms of a Markovian dynamics for the complex mode amplitudes of the expansion of the pressure field 
onto the guided modes of the unperturbed and non-dissipative waveguide.
This Markovian dynamics involves coupling terms between 
guided modes and mode-dependent dispersion and loss terms.
The coupling terms between guided modes come from the random perturbations in the water column.
The effective dispersion and loss terms come from two effects: the random perturbations in the water  column
induce coupling between the guided and radiative modes 
and the deterministic dissipation in the sediment layer induce an exponential decay of the guided mode amplitudes.
Both effects generate an irreversible, mode-dependent loss of energy carried by the guided modes.

The mathematical literature contains a lot of results on wave propagation in randomly perturbed waveguides 
motivated by underwater acoustics <cit.>.
Those results derive from first principles and make it possible to relate the coefficients of the effective Markovian model
to the physical parameters of the waveguide (in particular, the statistics of the fluctuations of the index 
of refraction and the complex acoustic impedance of the sea bottom).
The mathematical statements also make it clear in which sense the Markovian model approximates the random dynamics
of the mode amplitudes in the waveguide,
because there are subtle effects that follow from the fact that the results are first established in a weak topology, 
as explained in Appendix <ref>.
Comparisons between theory and numerical simulations show good agreement, 
both for direct and synthetic inverse problems <cit.>,
but there is not so far much comparison between the detailed mathematical predictions and real experiments.
On the other hand, the physical literature contains many theoretical results essentially based on coupled mode equations <cit.>.
However, the equations are derived ad hoc and it is not always straightforward to relate the coefficients
of the effective equations to the physical
and statistical parameters of the medium. 
This is problematic as we have in mind to use such results to solve an inverse problem.

As we will see in this paper, the effective attenuation of the mode amplitudes plays a key role.
In the physical literature, modal attenuation coefficients are introduced directly in the coupled mode equations without derivation from first principles
 <cit.>.
On the other hand most of mathematical studies do not take into account the attenuation of the bottom layer.
However, even though dissipation is weak, it plays an important role in long-range propagation.
For the direct problem point of view, the equipartition regime that is well characterized by Gaussian statistics <cit.>
and a scintillation index (relative intensity variance) 
close to one in the absence of dissipation looks very different in the presence of weak attenuation and it
may give rise to high fluctuations of the intensity. This was first pointed out by Creamer <cit.>.
For the inverse problem point of view, signals measured on an hydrophone array can be processed
to extract these attenuation coefficients, which in turn allow to get information on the medium,
in particular, the bottom properties, as we will see in Section <ref>.
We exploit data from at-sea experiments based on DGA's ALMA (Acoustic Laboratory for 
Marine Applications) system <cit.>.
We use data collected by a system consisting of four 32-element, vertical-line hydrophone array and a moored pinger located far away (at 9 km)
and  transmitting time-harmonic signals between 2 and 13 kHz with a three-minute repetition rate.
The measurements were carried out in November (7th-17th) 2016 near the shores of North East Corsica <cit.>.
The signals recorded by the hydrophones can be processed by cross correlation calculations
to estimate the properties of the medium. We explain this procedure  and report its results  in this paper.

The paper is organized as follows.
The direct problem is analyzed in Sections <ref>-<ref>. Section <ref> is a review of the modal decomposition
of the sound pressure in a homogeneous, non-dissipative waveguide. Section <ref> describes the Markovian dynamics
of the mode amplitudes in a random, dissipative waveguide.
The inverse problem is formulated and solved using experimental data in Section <ref>.




§ WAVE PROPAGATION IN WAVEGUIDES


Our model consists of a two-dimensional waveguide with range axis denoted by x ∈ and transverse coordinate
denoted by z ∈ [0,+∞). We suppose that the depth of the sea is constant and equals z_ b. When z ∈ [0,z_ b] the medium is water, when z> z_ b it becomes sediments. A point-like source at a fixed position (x,z) = (0,z_0) transmits a time-harmonic signal at frequency ω 
which is collected by a vertical array of receivers (hydrophones) at x=x_ a.

The Helmholtz equation for the acoustic pressure p̂(x,z) writes:

    [  
     ρ(x,z) ∇·( 1/ρ(x,z)∇) + ω^2/c(x,z)^2] 
    p̂(x,z) =  δ(x)δ(z-z_0),

for x ∈,   z ≥ 0,
where ∇=(∂_x,∂_z), c(x,z), resp. ρ(x,z), is the speed, resp. the density, at position (x,z).  
We consider that ρ is stepwise constant and equal to ρ_ w  in the water and 
ρ_ s in the sediments.
Therefore we have

    [  (∂_x^2 + ∂_z^2 ) +  ω^2/c(x,z)^2] 
    p̂(x,z)  = δ(x)δ(z-z_0),

for x ∈,   z ∈ (0,z_ b)∪(z_ b,+∞).


Remark.
The model (<ref>) can be derived from a more realistic three-dimensional situation, in which 
the sound pressure field P̂ satisfies in cylindrical coordinates 

    [ ∂_r^2  +1/r∂_r +1/r^2∂_θ^2  +∂_z^2  +
    ω^2/c(r,z)^2] 
    P̂ = 1/2π rδ(r)δ(z-z_0) .

The solution is radially symmetric and, 
if we neglect a near field factor of the form p̂/r^5/2 then the scaled pressure field 
p̂(r,z) = √(r)P̂(r,z)
satisfies (<ref>).


The differential equation (<ref>) is endowed with the following boundary conditions: 

    [ ∙; ∙ ]


We are interested in solutions of (<ref>) such that

    p̂ 1_(0,+∞)(x) 
    ∈ C^0( (0,+∞),H_0^1(0,+∞)∩ H^2(0,+∞))
    ∩ C^2( (0,+∞),L^2(0,+∞)) .


Let us now describe the results we obtain in this model. 
We compute in Section <ref> the pressure 
for an ideal (homogeneous) waveguide in an exact form. 
We analyze how it behaves when the sound speed in the water is no longer constant but perturbed 
by some small noise and the sediments are weakly dissipative. 
The wave modes interact with each other and an asymptotic regime in the small noise and large distance limit is studied in Section <ref>. 
In particular we compute correlations of the recorded pressure signals for hydrophones 
located in a vertical segment  for exponentially decaying correlations of the medium in Section <ref>.



§ HOMOGENEOUS WAVEGUIDE


In this section, we consider a wave speed c(x,z) which is constant both in the water and in the sediments:

    c_0(z) = c_ w 1_[0,z_ b](z) + c_ s 1_( z_ b,+∞) (z) ,
 
with c_ s > c_ w (Pekeris waveguide). There is no dissipation, no fluctuation along the x-axis.
The analysis of the perfect waveguide is classical <cit.>
but we include it here for the sake of completeness.
Let us  introduce the Helmholtz operator

    H = 
    ρ_0(z) ∂_z ρ_0(z)^-1∂_z + ω^2/c_0(z)^2  ,

where ρ_0(z) = ρ_ w 1_[0,z_ b](z) + ρ_ s 1_( z_ b,+∞) (z),
with boundary conditions:

    [ ∙; ∙ ]

We also denote k_0(z) = k_ w 1_[0,z_ b](z) + k_ s 1_( z_ b,+∞)(z),
with k_ w = ω/c_ w and k_ s = ω/c_ s.
The Helmholtz operator  is self-adjoint with respect to the scalar product defined on L^2(^+,ρ_0^-1) by:

    (ϕ_1,ϕ_2)    := ∫_0^∞ρ_0^-1(z) ϕ_1(z)ϕ_2(z)dz 
       = ρ_ w^-1∫_0^z_ bϕ_1(z)ϕ_2(z)dz +  ρ_ s^-1∫_z_ b^+∞ϕ_1(z)ϕ_2(z)dz.

The Helmholtz operator has a spectrum of the form

    (-∞,  k_ s^2) ∪{β_N^2,…,β_1^2}  ,

where the N modal wavenumbers β_j are positive and 
k_ s^2<β^2_N < ⋯ <β_1^2 <k_ w^2.

Discrete spectrum.
The jth eigenvector associated to the eigenvalue β_j^2 is

    ϕ_j(z) = 
    {[                A_j sin (σ_j z/z_ b)  0≤ z≤ z_ b; A_j sin (σ_j ) exp(-ζ_j (z-z_ b)/z_ b)  z≥ z_ b ].

where 

    σ_j = z_ b√(k_ w^2 -β_j^2),   ζ_j =z_ b√(β_j^2-k_ s^2)  ,

and

    A_j^2 = 2/z_ b/1/ρ_ w ( 1-sin (2σ_j)/2σ_j)
    +1/ρ_ ssin^2 (σ_j)/ζ_j  .

The σ_j's are the solutions in (0, √(k_ w^2-k_ s^2) z_ b) of

    tan (σ) = - σ/√( z_ b^2 (k_ w^2-k_ s)^2-σ^2)ρ_ s/ρ_ w ,

and we denote by N the number of solutions.


Continuous spectrum.
For γ∈ (-∞,k_ s^2), the improper eigenvector has the form:

    ϕ_γ(z) = 
    {[                                                                     A_γsin (η_γ z/z_ b)  0≤ z≤ z_ b; A_γ[ sin (η_γ ) cos(ξ_γ (z-z_ b)/z_ b)
    +ρ_ s/ρ_ wη_γ/ξ_γcos(η_γ) sin(ξ_γ (z-z_ b)/z_ b)
    ]
      z≥ z_ b ].

where

    η_γ=z_ b√(k_ w^2-γ),  ξ_γ=z_ b√(k_ s^2-γ)

and

    A_γ^2 = ξ_γρ_ s z_ b/π (ξ_γ^2 sin^2(η_γ) +
    ρ_ s^2/ρ_ w^2η_γ^2 cos^2(η_γ))  .

We remark that ϕ_γ does not belong to L^2(^+, ρ_0^-1), but 
( ϕ_γ,ϕ) can be defined for any test function ϕ∈ L^2(^+,ρ_0^-1) as

    ( ϕ_γ,ϕ) = lim_M → +∞∫_0^M ϕ_γ(z) ϕ(z) ρ_0(z)^-1 dz
    ,

where the limit holds on L^2(-∞,k_ s^2).

Completeness.
We have for any ϕ∈ L^2(^+, ρ_0^-1):

    ( ϕ,ϕ)  = ∑_j=1^N | ( ϕ_j ,ϕ) |^2
    +∫_-∞^k_ s^2| ( ϕ_γ,ϕ)|^2 dγ   .

The map which assigns to every element of L^2(^+, ρ_0^-1) the coefficients of its spectral decomposition

    ϕ↦( ( ϕ_j,ϕ)  ,j=1,…,N, 
     (ϕ_γ,ϕ), γ∈ (-∞,k_ s^2))

is an isometry from L^2(^+, ρ_0^-1) onto ^N × L^2(-∞,k_ s^2).
There exists a resolution of the identity Π such that for any ϕ∈ L^2(^+ , ρ_0^-1) and -∞≤ r ≤ r' ≤ +∞:

    Π(r,r')(ϕ) = ∑_j=1^N ( ϕ_j ,ϕ) ϕ_j  1_(r,r')(β_j^2)
    +
    ∫_r^min(k_ s^2,r')( ϕ_γ, ϕ) ϕ_γ dγ 1_(r,+∞)(k_ s^2) ,

and for any ϕ in the domain of H,

    Π(r,r')( Hϕ)  = ∑_j=1^N β_j^2 ( ϕ_j ,ϕ) ϕ_j  1_(r,r')(β_j^2)
    +
    ∫_r^min(k_ s^2,r')γ( ϕ_γ, ϕ) ϕ_γ dγ 1_(r,+∞)(k_ s^2) .


Modal decomposition.
Any solution of the Helmholtz equation in homogeneous medium can be expanded as 

    p̂(x,z) = ∑_j=1^N p̂_j(x) ϕ_j(z)
    +∫_0^k_ s^2p̂_γ (x) ϕ_γ(z) dγ
    +∫_-∞^0 p̂_γ (x) ϕ_γ(z) dγ .

The first modes are guided, the second ones  are radiating, the third ones are evanescent.
Indeed, the complex mode amplitudes satisfy

    ∂_x^2 p̂_j +β_j^2 p̂_j   =0,     j=1,…,N,
    ∂_x^2 p̂_γ +γp̂_γ   =0,   γ∈ (-∞,k_ s^2) .

Therefore, if the source is in the plane x=0, at (0,z_0), as in (<ref>), we have for x>0:

    p̂(x,z) =   ∑_j=1^N â_j,0/√(β_j) e^i β_j xϕ_j(z)
    +∫_0^k_ s^2â_γ,0/γ^1/4 e^i √(γ) xϕ_γ(z) dγ
       +∫_-∞^0  â_γ,0/|γ|^1/4   e^- √(|γ|) xϕ_γ(z) dγ
    ,

with the modal amplitudes determined by the source:

    â_j,0 =     √(β_j)/2ϕ_j(z_0)  ,   j=1,…,N,
    â_γ,0 =    γ^1/4/2ϕ_γ(z_0),  γ∈ (0,k_ s^2) ,
    â_γ,0 =    |γ|^1/4/2ϕ_γ(z_0) ,  γ∈ (-∞,0) .






§ RANDOM AND DISSIPATIVE WAVEGUIDE


Here we assume that the waveguide is weakly randomly perturbed and weakly dissipative:

    c^2(x,z) = c_0^2(z)/1+V(x,z) .

The perturbation V(x,z) has two components: random velocity perturbation in water
 and deterministic constant dissipation in sediments:

    V(x,z) =  ν(x,z)  1_(0,z_ b)(z)
    + i  ν_ s 1_ (z_ b,+∞)(z),

where ν(x,z) is a zero-mean random process that describes the relative fluctuations
of the propagation speed in water and ν_ s>0 models the damping in the sediments.



 §.§ The coupled mode equations

The solution of the perturbed Helmholtz equation

    [(∂_z^2 + ∂_x^2 )+ ω^2/c_0^2(z) + ω^2/c_0^2(z) V(x,z)] p̂( x,z) =  
    δ(x)δ(z-z_0) ,

for x∈, z∈ (0,z_ b)∪(z_ b,+∞), 
with the boundary conditions (<ref>),
can be expanded as (<ref>) and the complex mode amplitudes satisfy the coupled equations:

    ∂_x^2 p̂_j +β_j^2 p̂_j   = -ω^2 ∑_l=1^N C_jl(x) p̂_l 
    -ω^2 ∫_-∞^k_ s^2 C_jγ' (x) p̂_γ' dγ',     j=1,…,N,
    ∂_x^2 p̂_γ +γp̂_γ   = -ω^2 ∑_l=1^N C_γ l(x) p̂_l 
    -ω^2 ∫_-∞^k_ s^2 C_γγ' (x) p̂_γ' dγ',  
      γ∈ (-∞,k_ s^2)  ,

with

    C_jl(x) =   (ϕ_j,ϕ_l V(x,·) c_0^-2) ,
    
    C_jγ'(x) =   (ϕ_j,ϕ_γ' V(x,·) c_0^-2) ,
    
    C_γ l(x) =   (ϕ_γ,ϕ_l V(x,·)c_0^-2) ,
    
    C_γγ'(x) =   (ϕ_γ,ϕ_γ' V(x,·) c_0^-2) .

The coupling term C_jl has two components:

    C_jl(x) =  C^ w_jl(x) + i  ν_ s ρ_ s^-1 c_ s^-2∫_z_ b^+∞ϕ_j(z) ϕ_l(z) dz ,

where

    C^ w_jl(x) = c_ w^-2ρ_ w^-1∫_0^z_ bν(x,z) ϕ_j(z) ϕ_l(z) dz ,

and similarly for the other coupling terms.



We introduce the amplitudes of the generalized right and left-going mode amplitudes â_j and b̂_j
for the guided and radiating modes:

    p̂_j(x) =   1/√(β_j)( â_j (x) e^iβ_j x +b̂_j(x) e^- iβ_j x),    j=1,…,N,
    ∂_x  p̂_j(x) =    i√(β_j)( â_j (x) e^iβ_j x -b̂_j(x) e^- iβ_j x),   j=1,…,N,
    p̂_γ(x) =   1/γ^1/4( â_γ (x) e^i√(γ) x 
    +b̂_γ(x) e^- i√(γ)x),  γ∈ (0,k_ s^2) ,
    ∂_x  p̂_γ(x) =    iγ^1/4( â_γ (x) e^i√(γ) x 
    -b̂_γ(x) e^- i√(γ)x),  γ∈ (0,k_ s^2) .

They satisfy the following coupled differential equations:

    d/dxâ_j(x) =   i ω^2/2∑_l =1^N C_jl(x)/√(β_j β_l)(â_l( x) e^i(β_l-β_j) x + b̂_l( x) e^-i (β_l+β_j) x) 
       + i ω^2/2∫_0^k_ s^2C_j γ'(x)/√(β_j √(γ'))(â_γ'( x) e^i(√(γ')-β_j) x + b̂_γ'(x) e^-i (√(γ')+β_j) x) dγ' 
       +i ω^2/2∫_-∞^0C_j γ'(x)/√(β_j)p̂_γ'( x) e^- i β_j x dγ'  ,


    d/dxb̂_j( x) =    - i ω^2/2∑_l =1^N C_jl(x)/√(β_j β_l)(â_l( x) e^i(β_l+β_j) x + b̂_l(x) e^-i (β_l-β_j) x) 
       - i ω^2/2∫_0^k_ s^2C_j γ'(x)/√(β_j √(γ'))(â_γ'(x) e^i(√(γ')+β_j) x + b̂_γ'( x) e^- i (√(γ')-β_j) x) dγ' 
       -i ω^2/2∫_-∞^0C_j γ'(x)/√(β_j)p̂_γ'(x) e^i β_j x dγ'  .

The radiating  mode amplitudes â_γ (x) and b̂_γ (x), γ∈ (0,k_ s^2),
satisfy the same equations upon replacing j by γ in the formulas above.
The evanescent mode amplitudes p̂_γ(x), γ∈ (-∞,0), satisfy the differential equations

    ∂_x^2 p̂_γ +γp̂_γ = 
    -  g_γ(x) -ĝ_γ^ ev(x) ,    x  0,

where ĝ_γ(x) and  ĝ_γ^ ev(x)  
are given by

    ĝ_γ(x) =   ω^2
    ∑_l'=1^NC_γ l'(x)/√(β_l')[ a_l'(x) e^i β_l' x
    + b_l'(x) e^- iβ_l' x]
       +   ω^2 
    ∫_0^k_ s^2C_γγ'(x)/√(γ')[ a_γ'(z) e^i  √(γ') x
    + b_γ'(x) e^- i √(γ') z]dγ',  
    ĝ^ ev_γ(x) =   ω^2  ∫_-∞^0 C_γγ'(x) p̂_γ'  (x)
    dγ' .


We wish to understand the evolution of the complex mode amplitudes in the regime where the perturbations are weak and 
the propagation distance is large. We assume 

    V(x,z) =  ν(x,z)  1_z∈ (0,z_ b)
    + i ^2 ν_ s 1_z∈ (z_ b,+∞),

where  is a small dimensionless parameter that expresses the fact that the random perturbations are weak and dissipation is even weaker.
We suppose that we have a perfect (homogeneous and non-dissipative) waveguide outside the region x∈(0,x_∞/^2) for some x_∞>0.
We carry out an asymptotic study in the small  limit. 
In this framework the evanescent mode amplitudes can be expressed in terms of the propagating mode amplitudes <cit.>
and the differential equation for â_j^(ω,x):=â_j(ω,x/^2) can be written as:

    d/dxâ_j^(x)    = i ω^2/2∑_l =1^N ^-1 C^ w_jl(x/^2) + i ν_ sρ_ s^-1 c_ s^-2∫_z_ b^∞ϕ_j(z) ϕ_l(z) dz/√(β_j β_l)
           ×  (â_l^(x) e^i(β_l-β_j) x/^2+   b̂_l^(x) e^-i (β_l+β_j) x/^2) + O(1).

In the asymptotic framework “weak fluctuations, weak dissipation, large propagation distance", 
we carry out a multiscale analysis similar as the one reported in Ref.  gomez.
We find that the random fluctuations of the propagation speed in water give the standard coupling terms
between the right-going mode amplitudes of the guided modes. They also induce loss towards the radiating modes.
Moreover, the dissipation in the sediments gives rise to additional damping terms in the evolution equations
of the complex mode amplitudes.
The coupled system of equations for the complex mode amplitudes is of the form d a^/dx = ^-1 F(x/^2,  a^) +  G(x/^2,  a^). 
The small  limit of such differential equations are studied in the book  FGPSbook (Chapters 6 and 20). 
The full result for the complex mode amplitudes is expressed in Appendix <ref>.
From this result we get that the process (|â_j^(x)|^2)_j=1^N converges towards a Markov process P(x)=(P_j(x))_j=1^N whose infinitesimal generator ℒ_ P writes:

    ℒ_ P = ∑_j ≠ lΓ_lj[P_l P_j (∂/∂ P_j - ∂/∂ P_l) ∂/∂ P_j + (P_l - P_j) ∂/∂ P_j] -  ∑_j=1^N Λ_j P_j ∂/∂ P_j ,

with

    Λ_j =   ∫_0^k_ s^2ω^4/2√(γ)β_j∫_0^∞ [ C^ w_jγ(0)C^ w_jγ(x) ] cos((√(γ)-β_j)x) dx dγ
       +
    ν_ sω^2/β_j c_ s^2 ρ_ s∫_z_ b^∞ϕ_j(z)^2 dz  , 
    Γ_lj=   ω^4/2 β_j β_l∫_0^∞ [C^ w_jl(0)C^ w_jl(x)] cos((β_l-β_j)x) dx    j ≠ l,

where

    [ C^ w_jγ(0)C^ w_jγ(x) ]
    =   ρ_ w^-2  c_ w^-4∫_0^z_ b∫_0^z_ bϕ_j ϕ_γ(z)  [ν(0,z)ν(x,z')] ϕ_j ϕ_γ(z') dz dz' ,
     [C^ w_jl(0)C^ w_jl(x)]  =   ρ_ w^-2  c_ w^-4∫_0^z_ b∫_0^z_ bϕ_j ϕ_l(z)  [ν(0,z)ν(x,z')] ϕ_j ϕ_l(z') dz dz' ,

and we have ∫_z_ b^∞ϕ_j(z)^2 dz = A_j^2 sin^2(σ_j)  z_ b/(2ζ_j).
This means that, for any test function F:ℝ^N →ℝ,
the expectation [ F(  P(x) )] satisfies the Kolmogorov forward equation

    ∂_x  [ F(  P(x) )] =  [ L_ PF(  P(x) )] .

Equivalently, the probability density function p(x, P) of the random vector P(x) satisfies the Kolmogorov backward
(or Fokker-Planck) equation
∂_x p=   L_ P^* p ,
where L_ P^* is the adjoint of L_ P.


From the form of the generator L_ P, one can establish that the nth-order moments 
of the mode powers satisfy closed equations.
In particular, using (<ref>) and (<ref>) we find that the mean mode powers

    Q_j(x)=
     [P_j(x)]  = lim_→ 0[ |â_j^(x)|^2]

satisfy the closed system of equations

    ∂_x Q_j = -Λ_j Q_j +∑_l=1^N Γ_lj(Q_l-Q_j)  ,
      Q_j(0)=|â_j,0|^2 .

The form of these coupled-mode equations is well-known <cit.> although the mode-dependent attenuation term Λ_j 
was usually introduced heuristically so far.
The solution explicitly writes:

    Q(x) = exp( A x)  Q(0)  ,

with the matrix A defined by (δ_jl is the Kronecker symbol and Γ_jj = -∑_l'≠ jΓ_jl'):

    A := (Γ_jl - Λ_j δ_jl)_j,l=1^N .



 


 §.§ Computation of the coefficients of the matrix A

The goal of this paragraph is to show how to get closed-form expressions of the coefficients of
 the matrix A when the correlation of the perturbation decreases exponentially 
as a function of the horizontal distance between the two points, i.e.

    [ν(x,z) ν(x',z')] = σ^2 exp(-|x-x'|/ℓ_ h) R(z,z') ,

where ℓ_ h is the horizontal correlation radius of the random fluctuations of the index of refraction.
We can compute Γ_jl for j ≠ l from:

    ∫_0^∞ [C^ w_jl(0)C^ w_jl(x)] cos((β_l-β_j)x) dx 
          = σ^2/ρ_ w^2 c_ w^4ℓ_ h/1+ (β_l - β_j)^2 ℓ_ h^2∫_0^z_ b∫_0^z_ b R(z,z') ϕ_j ϕ_l(z) ϕ_j ϕ_l(z') dz dz'

and similarly we can compute Λ_j from the same quantities upon substitution γ for l.
By (<ref>) we have

    Γ_lj=   ω^4/2 β_j β_lσ^2/ρ_ w^2 c_ w^4ℓ_ h/1+ (β_l - β_j)^2 ℓ_ h^2∫_0^z_ b∫_0^z_ b R(z,z') ϕ_j ϕ_l(z) ϕ_j ϕ_l(z') dz dz' .

Replacing ϕ_j for all j by their expression, and using the notation k_ wj := √(k_ w^2 -β_j^2), we obtain: 

    ∫_0^z_ b∫_0^z_ b R(z,z') ϕ_j ϕ_l(z) ϕ_j ϕ_l(z') dz dz' 
     = A_j^2 A_l^2/4(S(k_ wj - k_ wl, k_ wj - k_ wl) 
        + S(k_ wj + k_ wl, k_ wj + k_ wl)  - S(k_ wj - k_ wl, k_ wj + k_ wl) - S(k_ wj + k_ wl, k_ wj - k_ wl)),

where 

    S(k, k') :=  ∫_0^z_ b∫_0^z_ b R(z,z') cos(k  z) cos(k'  z') dz dz'.

When R has the exponential form R(z,z') = exp(-|z-z'|/ℓ_ v)/2, then
the variance of the fluctuations is σ^2/2 and S(k,k')
has a closed-form expression. 
We will use this particular form of the correlation function of the medium in Section <ref> because it gives simple expressions
for the coefficients Γ_jl and Λ_j, which is convenient for the resolution of the inverse problem
which requires many evaluations of such coefficients for different values σ, ℓ_ v,ℓ_ h of the statistics of the random medium 
and ρ_ s, ν_ s of the sea bottom.
We could as well use the  Garrett-Munk correlation function at the expense of some computational overburden <cit.>,
but the impact of the exact form of the correlation function turns out to be weak.




 §.§ Pressure field correlations


We now compute the correlation of the received signal on a vertical segment at a fixed horizontal distance x_ a. When x_ a≫ 1, we can use the asymptotic study from above and the decomposition of p̂( x_ a, z) in the basis of eigenvectors ϕ_j is: 

    p̂( x_ a, z) ≃∑_j=1^N 1/√(β_j)â_j (x_ a)  ϕ_j(z).

Therefore

    [p̂( x_ a, z)p̂( x_ a, z') ] 
       = ∑_j=1^N  [|â_j(x_ a)|^2]/β_jϕ_j(z)ϕ_j(z')+∑_j≠ l [â_j(x_ a)â_l(x_ a) ]/√(β_j β_l)ϕ_j(z)ϕ_l(z').

The cross second moments [â_j â_l  ] for j≠ l decay exponentially with the propagation distance x
as shown in Ref.  GP07 (see also Ref.  colosi12) and they can be neglected as soon as the propagation distance x is larger than the scattering mean free path,
that depends on the coefficients Γ_jl. This is roughly speaking the distance beyond which the phase of the 
field has a random part which variance of the order of or larger than one, so that the coherent field is vanishing <cit.>.
The second moments [|â_j |^2 ] were computed
 in the previous paragraph using the exponential of the matrix A:

    [|â_j(x_ a)|^2] =  ∑_l=1^N (exp( A x_ a))_jl |â_l,0|^2,

where â_l,0 is given by (<ref>).
We can now write the correlation between the sound pressure signals recorded by two receivers which are separated by y along the vertical segment at distance x_ a and of depth [z_m,z_M] where 0 <z_m < z_M < z_ b.
For all y ∈ [0,z_M -z_m], the spatial correlation C_x_ a(y) writes:

    C_x_ a(y)    := 1/z_M-z_m -y∫_z_m^z_M-y[p̂( x_ a, z)p̂( x_ a, z+y) ] dz 
       = 1/z_M-z_m -y∑_j  [|â_j(x_ a)|^2]/β_jA_j^2/2(cos(k_ wj y) (z_M-z_m-y) 
          - sin(k_ wj(2 z_M-y)) - sin(k_ wj(2 z_m+y))/2 k_ wj).






 §.§ Equipartition regime

By (<ref>) the mean mode powers satisfy

    Q_j(x) x → +∞≃
    c_V V_j  exp( -λ  x ) (1+o(1))  ,

where ( V ,-λ ) is the first eigenvector/eigenvalue of the matrix
A=Γ -Φ,
with Γ_jl given by (<ref>) for j≠ l,  Γ_jj = -∑_l'≠ jΓ_jl', 
Φ_jl = Λ_j δ_jl,
and

    c_V = ∑_j=1^N V_j  |â_j,0|^2  .

Note that:


  * The coefficients of V have all the same sign (so we can assume that they are nonnegative).

Indeed, the solution of the equation ∂_x  Q = (Γ-Φ)  Q starting from any
Q(x=0) with nonnegative coefficients must have nonnegative coefficients since the exponential of a matrix whose off diagonal entries are nonnegative has nonnegative coefficients.
As Q(x) is equivalent to ( V^T Q(x=0)) 
exp(-λ  x)  V, the coefficients of V must have the same sign. We deduce as well that c_V has the same sign as V.

  * λ≥0.

Indeed, we have (Γ -Φ )   V = - λ V. By projecting onto
V^(0) = (1/√(N))_j=1^N, which is in the kernel of Γ, we get
V^(0)^T Φ V = λ V^(0)^T  V.
Since Φ is diagonal with nonnegative coefficients and the coefficients of V have the same sign and
cannot be all zero, we find λ  = ∑_j Λ_j V_j / ∑_j V_j ≥ 0.





In the following we discuss cases with zero or weak dissipation, where explicit expressions can be obtained.
Remember that the effective dissipation is the sum of the two effects: power leakage from the guided modes to the radiating modes and
dissipation in the sediments.

No effective dissipation.
If there is no effective dissipation Φ = 0, 
then the first eigenvector/eigenvalue ( V^(0),-λ^(0)) of the matrix
Γ is 

    V^(0) = (1/√(N))_j=1^N,   λ^(0)=0 ,

which gives the standard equipartition result <cit.>:

    Q_j(x) x → +∞⟶1/N∑_j=1^N  |â_j,0|^2   .


Weak effective dissipation.
We next consider the case when the effective dissipation is weak, 
that is to say,  the matrix Φ is much smaller
  than the matrix Γ, with a typical ratio of the order of δ.
We then assume that
Λ_j = δΛ_j^(1)  ,
with δ≪ 1.
Then we can write Φ  =δΦ^(1) 
with Φ_jl^(1) =
Λ_j^(1)δ_jl
and Γ  = Γ^(0)
and the first eigenvector/eigenvalue ( V ,-λ)  of the matrix
Γ -Φ can be expanded
as

    V  =  V^(0) +δ V^(1) +O(δ^2)
      ,  λ= δλ^(1) +δ^2 λ^(2) +O(δ^3),

with

    λ^(1)   =    V^(0)^T Φ^(1) V^(0)
    =
    1/N∑_j=1^N Λ_j^(1),
    λ^(2)   =    V^(0)^T Γ^(0) V^(1) ,

and V^(1) is solution of Γ^(0) V^(1) = (Φ^(1) - λ^(1) )
 V^(0) and is orthogonal to V^(0).
If, for instance, Γ_jl≡Γ >0 for all j≠ l, then

    V^(1) = -1/Γ N^3/2(Λ_j^(1))_j=1^N

and

    λ^(2) = 
    - 1/Γ N^2∑_j=1^N (Λ_j^(1) -λ^(1))^2 .





 §.§ Fluctuation analysis

By (<ref>) we find that the second-order moments of the mode powers 

    R_jl  (x)= [P_j(x)P_l(x)] ,   j,l=1,…,N,

satisfy the closed equations

    ∂_x R_jj = - 2 Λ_j R_jj  +∑_n ≠ jΓ_jn (4R_jn -2R_jj )  ,
       ∂_x R_jl  = - (2Γ_jl+ Λ_j +Λ_l) R_jl  +
    ∑_n ≠ lΓ_ln (R_jn -R_jl )  + ∑_n ≠ jΓ_jn (R_nl -R_jl )  ,

for j≠ l. 
This system has the same form as the one found in the literature dedicated to coupled mode theory <cit.>.
The initial conditions are R_jl (0)=|â_j,0|^2|â_l,0|^2.
Let us introduce S = (S_jl)_1≤ j≤ l≤ N defined by

    S_jl = 
    {[ R_jl+R_lj  j  <  l,;      R_jj    j=l  . ].

The S_jl's satisfy the system

    ∂_x S_jl =    - (Λ_j+Λ_l)S_jl
    +∑_n ∉{j,l}[ Γ_ln (S_jn-S_jl)+Γ_jn(S_nl-S_jl)] 
       
    +2Γ_jl 1_j ≠ l (S_jj+S_ll -2S_jl) ,

with the convention that whenever S_jl occurs with j>l, it is replaced by S_lj.
This can be written in the form ∂_x  S  = (Θ 
- Ψ )  S.
The linear operator Ψ (that depends on the Λ_j's) is diagonal and the linear operator Θ
(that depends on the Γ_jl's)  is self-adjoint:  for any T and T, we have

    ∑_j≤ l  (Θ T)_jlT_jl   = -∑_j≤ lΘ_jl,jl T_jlT_jl
    +
    ∑_j <l , n ∉{j,l}( Γ_ln T_jnT_jl +Γ_jnT_nlT_jl) 
          +∑_j ≠ n(Γ_jn T_jnT_jj +Γ_jn T_njT_jj) 
    +2 ∑_j < l( Γ_jl T_jjT_jl + Γ_jl T_jjT_jl)
       =∑_j≤ l
    T_jl (Θ T)_jl 
    ,

because 2∑_j<l = ∑_j ≠ l.
As a consequence, Θ 
- Ψ can be diagonalized and 
we find that

    S(x) x → +∞≃
    c_W  Wexp( -μ x ) (1+o(1))  ,

where c_W is the projection on the first eigenvector W in the basis of eigenvectors of Θ-Ψ,
and

    c_W = 
    2 ∑_j< l W_jl  |â_j,0|^2|â_l,0|^2  +  ∑_j=1^N W_jj  |â_j,0|^4
    =
     ∑_j,l=1^N W_jl  |â_j,0|^2|â_l,0|^2  ,

with the convention that whenever W_jl occurs with j>l, it is replaced by W_lj.

No effective dissipation.
If there is no effective dissipation, then the first eigenvector/eigenvalue ( W^(0),-μ^(0)) of the matrix
Θ is 

    W^(0) = ( c_N )_1≤ j ≤ l ≤ N,   μ^(0)=0 ,

with c_N=√(2)/√(N(N+1)) 
and 

    S(x) x → +∞→
    c_W  W^(0)  .

As ∑_j≤ l S_jl(x) =∑_j,l R_jl(x) = (∑_j=1^N |â_j,0|^2)^2, we deduce 

    S(x) x → +∞→(∑_j=1^N |â_j,0|^2)^22/N(N+1) ,

and

    R(x) x → +∞→(∑_j=1^N |â_j,0|^2)^21+δ_jl/N(N+1).


Weak effective dissipation.
We next consider the case when the effective dissipation is weak, 
say Λ_j = δΛ_j^(1).
Then we can write Ψ =δΨ^(1) and Θ = Θ^(0)
and the first eigenvector/eigenvalue ( W,-μ)  of the matrix
Θ-Ψ can be expanded
as

    W =  W^(0) +δ W^(1) +O(δ^2)
      ,  μ= δμ^(1) +δ^2 μ^(2) +O(δ^3),

with

    μ^(1)   =    W^(0)^T Ψ^(1) W^(0)
    =
    c_N^2 ∑_j <l (Λ_j^(1)+Λ_l^(1))
    +
    c_N^2  ∑_j 2 Λ_j^(1)
       =   2/N∑_j=1^N  Λ_j^(1) = 2λ^(1),
    μ^(2)   =    W^(1)^T Θ^(0) W^(1) ,

and W^(1) is solution of Θ^(0) W^(1) =( Ψ^(1) -μ^(1))
 W^(0) and is orthogonal to W^(0).
If, for instance, Γ_jl≡Γ >0 for all j≠ l, then

    W_jl^(1)= - c_N/Γ N( Λ_j^(1)+ Λ_l^(1) -2λ^(1)) ,    j≤ l,

and

    μ^(2)  =
    ∑_j≤ l W_jl^(1) (  Θ^(0) W^(1) )_jl= - 2(N+2)/N^2(N+1) Γ∑_j=1^N (Λ_j^(1)-λ^(1))^2  .

Note that

    μ - 2λ   =   δ^2(μ^(2)  -2 λ^(2)) +O(δ^3)
    
       =   
    -2 δ^2/N^2(N+1) Γ∑_j=1^N (Λ_j^(1)-λ^(1))^2 
    +O(δ^3) .

is negative-valued.




Exponential growth of the intensity fluctuations.
It is a general feature that, for any matrix Γ and effective dissipation coefficients
Λ_j, we have
μ - 2λ≤ 0 (we have equality when there is no effective dissipation and this is a consequence of the forthcoming result (<ref>) and 
Cauchy-Schwarz inequality when there is dissipation).
The first two moments of the pointwise intensity 
|p̂(x,z)|^2 for large x are

    [ |p̂(x,z)|^2] =   ∑_j=1^Nϕ_j(z)^2/β_j c_V V_j e^- λ x ,
     [ |p̂(x,z)|^4 ] =   ∑_j,l=1^N ϕ_j(z)^2ϕ_l(z)^2/β_j β_l c_W W_jl e^- μ x .

Without dissipation we have the following result for the fluctuations of the pointwise intensity:

    [|p̂(x,z)|^4]/[|p̂(x,z)|^2]^2x →∞⟶2N/N+1  ,

which is equal to 2 when N ≫ 1,
and with dissipation

    [|p̂(x,z)|^4]/[|p̂(x,z)|^2]^2x →∞≃2N/N+1exp( - (μ - 2λ) x )  ,

that grows exponentially with the propagation distance (for very long distances, however, as μ - 2λ
is very small as shown above).
Eq. (<ref>) gives the expression of the exponential growth rate 
when dissipation is weak. In this regime the growth rate increases when the effective modal dissipation coefficients
become different from each other and decreases when the number of modes increases.






§ INVERSE PROBLEM


The goal of this section is to show that it is possible to estimate the statistical properties of the index of refraction
in the water column and the sea bottom properties from the sound pressure recorded by a vertical hydrophone array
and transmitted by distant time-harmonic sources.
This inverse problem can be formulated as a minimization problem that tries to match empirical quantities 
with theoretical ones that depend on the parameters to be estimated.
The theoretical model of the previous section shows that the correlation function of the sound pressure
has a non-trivial behavior that makes it possible to identify many relevant parameters as we explain below.
We first present the ALMA experimental data, then formulate the inverse problem, and finally estimate the model parameters. 




 §.§ ALMA 2016 experiment

ALMA is a series of at-sea experiments carried out by French DGA Naval Systems <cit.>. In 2016, 
the experiment took place in november 7th-17th near the shores of North East Corsica <cit.> (cf. Figure <ref>). 
There was a pinger immerged at fifty meters, and nine kilometers away, a passive array with 128 hydrophones (4 verticals arms of 32 hydrophones, the arms are 0.5 m apart from each other, the hydrophones are 0.15 m apart from each other) immerged at sixty meters. The sea bottom is relatively flat on the line between the pinger and the array (between 100 m and 115 m), and composed of sand and gravelly-sand. Sea was calm, with roughness height around 0.1 m (sea state 1).
The pinger transmitted, for about 2 hours and a half, a sequence consisting of several time-harmonic waves, with a repetition rate of three minutes.
The sequence was a train made of two-second time-harmonic waves at the following K=6 frequencies: 2 kHz, 5 kHz, 7 kHz, 9 kHz, 11 kHz, 13 kHz. Acquisition was sampled at 48 kHz. The Fourier coefficients for each frequency were extracted using a Hann window function with a duration of 1 s.




Each arm has the form of a vertical array at a distance x_ a of the source immerged at z_0. The hydrophones are regularly spaced between the depths z_m and z_M, with 0 < z_m < z_M < z_ b.
At frequency ω, for y ∈[ 0 , z_M - z_m ], we denote the spatial correlation between hydrophones at a distance y from each other by C_x_ a(y) (cf. Eq. (<ref>) for the theoretical expression predicted by our model). 
Then we define the correlation radius r as the 
half-width at half-maximum (i.e. the distance between hydrophones for which   the correlation is 1/2):

    C_x_ a(r) = 1/2 C_x_ a(0)  =1/2 .

From the experimental data at frequencies f_1, … , f_K (with ω=2π f) we can extract the experimental correlation radii r_ e(f_1), … , r_ e(f_K). For Φ a set of parameters of our theoretical model (assuming that we know the sound speed c_ w and density ρ_ w in water and the depth of the sea bottom),

    Φ = ( c_ s , ρ_ s , ν_ s , σ , ℓ_ v , ℓ_ h)   ,

we can compute the theoretical correlation radii r_ t(f_1,Φ) , … , r_ t(f_K,Φ). Then we can define the misfit function

    E(Φ) = ∑_k=1^K ( r_ t(f_k , Φ) - r_ e(f_k) )^2  ,

and we can determine the parameters Φ by minimizing over Φ the misfit function E(Φ).



 §.§ Estimation of the model parameters

We assume that sound speed is constant in water, taking value c_ w=1523 m/s from CTD (conductivity, temperature, and depth) measurements, and we set density of water ρ_ w at 1000 kg/m^3 (we could have taken a more precise value, but density has little influence on the model).  Here, we have x_ a = 9000 m. Both source and array are in the middle of the water column.
Figure <ref> presents the correlation radii obtained from experimental data for each frequency, and the theoretical correlation radii 
predicted by our model with the parameters determined by minimization of the misfit function  (<ref>):
    c_ s = 1630 m/s;
    α = 1.09 dB / wavelength;
    σ = 0.002 (relative fluctuation);
    ℓ_ v = 30 m;
    ℓ_ h = 100 m;
    ρ_ s = 1700 kg/m^3.
Parameter ν_ s is computed from the attenuation coefficient α expressed in decibel by wavelength. 
The values of the parameters seem compatible with experimental measurements carried out in similar environments <cit.>.
Figure <ref>  shows
 theoretical and experimental correlation functions at the different frequencies.





We carry out a basic (one-at-a-time) sensitivity analysis: we modify one by one the optimal values we have found, letting the others unchanged.
Results are in Figure <ref>. We see that c_ s has a noticable influence, as a 20 m/s change in sediment sound speed significantly modifies the  correlation radii. σ is also important, as dividing or multiplying it by two has a strong influence. α and ℓ_ v have less influence, but are still significant. Finally, for reasonable physical changes of ρ_ s and ℓ_ h, we do not see noticeable effects in the results, whatever the frequency is.


The inspection of the numerical values of the two terms that determine Λ_j (cf. Eq. (<ref>)) give interesting information.
We can see that for each frequency, the radiative terms are approximately ten times smaller than the attenuation ones. 
So we can claim that, in the experimental configuration addressed in this paper, radiation effects can be neglected.



Finally, if we denote by  (z_n)_1⩽ n ⩽ N_ h the depths of the N_ h hydrophones of the array, 
we can compare the theoretical frequency-dependent scintillation indices defined by

    S_ t = 1/N_ h∑_n=1^N_ h[|p̂(x_ a,z_n)|^4] - [|p̂ (x_ a,z_n)|^2]^2/[|p̂(x_ a,z_n)|]^2 ,

with the experimental values S_ e determined by empirical averages instead of expectations.
Table <ref> shows experimental scintillation indices compared to theoretical ones, with the optimal parameters Φ 
determined from the correlation radii. In the experimental and theoretical configurations scintillation indices are approximately equal to one. One should 
have longer propagation distances to observe scintillation indices larger than one, as already suggested in Ref.  creamer.






§ CONCLUSION

This paper proposes a complete description of the statistics of the mode amplitudes of the sound pressure in a shallow-water waveguide.
The effective parameters (frequency- and mode-dependent attenuation, dispersion, and coupling) are identified from first principles
and expressed in terms of the statistical properties of the index of refraction of the water column and the sea bottom properties.
This theoretical analysis makes it possible to formulate an inverse problem for the estimation of the model parameters 
from the sound pressure recorded by a vertical  hydrophone array and transmitted by distant time-harmonic sources.
This inverse problem is solved using data collected during the ALMA 2016 experiment.

In the experimental configuration addressed in this paper, it turns out that radiation effects can be neglected (compared to dissipation in the sediments) 
and non-Gaussian scintillation effects can be neglected as well.
The extraction of the  frequency-dependent correlation radii of the recorded sound pressure signals 
makes it possible to estimate different acoustic and geoacoustic parameters:
The sediment sound speed c_ s and the standard deviation of the index of refraction in the water column σ can be robustly estimated 
because they have strong effects on the correlation radii of the sound pressure. The attenuation parameter in the sediments α and the vertical  correlation radius of the index of refraction in the water column ℓ_ v can be estimated, because they have less influence, but are still significant. 
The sediment density ρ_ s and the horizontal correlation radius ℓ_ h are difficult to estimate because they 
have small effects of the correlation radii of the sound pressure.



§ ACKNOWLEDGEMENTS

We thank Dominique Fattaccioli and Gaultier Real from DGA Naval Systems
for interesting and stimulating discussions and for making us available the ALMA 2016 data.




§ THE EFFECTIVE MARKOVIAN SYSTEM FOR THE COMPLEX MODE AMPLITUDES


We assume that the reduced wavenumbers β_j  are distinct, for 1 ≤ j ≤ N. Then, for any x_∞>0, the  process 
( (â_j^(x) )_j=1^N, (â_γ^(x) )_γ∈ (0,k_ s^2))
converges in distribution in C^0([0,x_∞], ℂ^N× L^2(0,k_ s^2) ), where ℂ^N× L^2(0,k_ s^2)
is equipped with the weak topology,  to the Markov process  
( (a_j(x) )_j=1^N , (a_γ (x))_γ∈ (0,k_ s^2))
with infinitesimal generator
L=  L^1+ L^2+ L^3 ,
where L^j, 1 ≤ j ≤ 3, are the differential operators:

    L^1 =    1/2∑_j,l=1^NΓ_jl( 
    a_ja_j∂_a_l∂_a_l
    +
    a_la_l∂_a_j∂_a_j
    -
    a_ja_l∂_a_j∂_a_l
    -
    a_ja_l∂_a_j∂_a_l)
     1_ j≠ l
       +
     1/2∑_j,l=1^NΓ^1_j l( 
    a_ja_l∂_a_j∂_a_l
    +
    a_ja_l∂_a_j∂_a_l
    a_ja_l∂_a_j∂_a_l
    -
    a_ja_l∂_a_j∂_a_l)
    
    
       + 1/2∑_j=1^N( Γ_jj - Γ^1_jj)
    ( a_j∂_a_j + a_j∂_a_j)
    
    +i/2∑_j=1^NΓ^s_jj( a_j∂_a_j - a_j∂_a_j)  , 
     L^2
    =   
    -1/2∑_j=1^N ( Λ_j  +iΛ^s_j ) a_j∂_a_j 
    +
    ( Λ_j  - iΛ^s_j )  a_j∂_a_j  ,
     
     L^3
    =   
    i  
    ∑_j=1^Nκ_j( a_j∂_a_j 
    - a_j∂_a_j) .

In these definitions we use the  classical complex derivative:
if ζ=ζ_r+iζ_i, then ∂_ζ=(1/2)(∂_ζ_r-i ∂_ζ_i) and
∂_ζ =(1/2)(∂_ζ_r +i ∂_ζ_i),
and the coefficients of the operators (<ref>-<ref>) are defined for 
indices j ,l= 1, …, N, as follows: 

 For all j ≠ l, Γ_jl is defined by (<ref>) and 

    Γ^s_j l =   ω^4/2 β_jβ_l∫_0^∞[ C_jl^ w(0) C_jl^ w(x) ] sin[(β_l-β_j)x] dx .

For all j,l:

    Γ^1_jl =   ω^4/2 β_jβ_l∫_0^∞[ C_jj^ w(0) C_ll^ w(x) ]   dx.

For all j, Λ_j is defined by (<ref>) and

    Γ_jj =    -∑_l =1,l≠ j^NΓ_jl ,    Γ^s_jj = -∑_l =1,l≠ j^NΓ^s_jl  ,
    Λ_j^s =   ∫_0^k_ s^2ω^4/2√(γ)β_j∫_0^∞[ C_j γ^ w(0) C_jγ^ w(x) ] 
    sin[ (√(γ)-β_j)x] dx dγ, 
    κ_j =    ∫_-∞^0  ω^4/2√(|γ|)β_j∫_0^∞[ C_j γ^ w(0) C_j γ^ w(x) ] cos ( β_j x) e^-√(|γ|)xdx dγ.
 

Remarks: 

1) The convergence result holds in the weak topology. This means that we can only compute quantities of the 
form [ F( a_1,…,a_N, ∫_0^k_ s^2α_γ a_γ dγ)]
for any test function α∈ L^2(0,k_ s^2) and F:ℝ^N+1→ℝ,
which are the limits of [ F( â^_1,…,â^_N, ∫_0^k_ s^2α_γâ^_γ dγ)] as → 0. 

2) The generator L does not involve ∂_a_γ. Therefore 
(â_j^ (x))_j=1^N
converges in distribution in C^0([0,x_∞], ℂ^N  ) to
the Markov process  (a_j(x) )_j=1^N
with  generator L. The weak and strong topologies are the same  in ℂ^N,
so we can compute any moment of the 
form [ F( a_1,…,a_N)], which are the limits of [ F( â^_1,…,â^_N)].

3) L_1 is the contribution of the coupling between guided modes, 
which gives rise to power exchange between the guided modes;
L_2  is the contribution of the coupling between guided and radiating modes, 
which gives rise to power leakage from the guided modes to the radiating ones (effective diffusion)
and addition of frequency-dependent phases on the guided modes (effective dispersion);
L_2 also contains the effective mode-dependent term due to dissipation in the sediments;
L_3  is the contribution of the coupling between guided and evanescent modes,
which gives rise to additional phase terms on the guided modes (effective dispersion <cit.>).

4) If the generator L is applied to a test function that depends only on the mode powers (P_j = |a_j|^2 )_j=1^N,
then the result is a function that depends only on (P_j )_j=1^N. Thus, 
the mode powers (P_j(x) )_j=1^N define  a Markov process, with infinitesimal generator defined by (<ref>).

5) The radiation mode amplitudes  remain constant on L^2(0,k_ s^2), equipped with the weak topology, as → 0. However, this does not
describe the power ∫_0^k_ s^2 |â^_γ|^2 dγ transported by the radiation modes,  because the convergence does not hold in the strong topology of L^2(0,k_ s^2) so we do not have ∫_0^k_ s^2 |â^_γ|^2 dγ→∫_0^k_ s^2 |a_γ|^2 dγ as → 0.


99

beran
M. J. Beran and S. Frankenthal,
Volume scattering in a shallow channel,
J. Acoust. Soc. Am. 91, 3203-3211 (1992).

borcea15
L. Borcea, J. Garnier, and C. Tsogka, A quantitative study of source imaging in random waveguides, 
Commun. Math. Sci. 13, 749-776 (2015).

colosi09
J. A. Colosi and A. Morozov,  
Statistics of normal mode amplitudes in an ocean with random sound speed perturbations: 
Cross mode coherence and mean intensity,
J. Acoust. Soc. Am. 126, 1026-1035 (2009).

colosi12
J. A. Colosi, T. F. Duda, and A. K. Morozov,
Statistics of low-frequency normal-mode amplitudes in an ocean with random sound-speed perturbations: 
Shallow-water environments,
J. Acoust. Soc. Am. 131, 1749-1761 (2012).

creamer
D. Creamer, 
Scintillating shallow water waveguides,
J. Acoust. Soc. Am. 99, 2825-2838 (1996).

dozier
L. B.  Dozier and F. D.  Tappert, 
Statistics of normal-mode am- plitudes in a random ocean. I. Theory, 
J. Acoust. Soc. Am. 63, 353-365 (1978). 

dozier2
L. B.  Dozier and F. D.  Tappert, 
Statistics of normal-mode am- plitudes in a random ocean. II. Computations,
J. Acoust. Soc. Am. 64, 533-547 (1978).

dozier83
L. B. Dozier, 
A coupled mode model for spatial coherence of bottom-interacting energy,
in Proceedings of the Stochastic Modeling Workshop, 
edited by C. W. Spofford and J. M. Haynes, ARL-University of Texas, Austin TX, 1983.

alma15
D. Fattaccioli,
ALMA: a new experimental acoustic system to explore coastal and shallow waters, 
UACE 2015 Proceedings (Chania, Greece).


alma17c
D. Fattaccioli and G. Real,
The DGA "ALMA" Project: an overview of the recent improvements of the system capabilities and of the at-sea campaign ALMA-2016,
UACE 2017 Proceedings (Skiathos, Greece),
available at www.uaconferences.org/docs/2017_papers/649_UACE2017.pdf

flatte
S. M.  Flatté, R. Dashen,  W. H. Munk, K. M. Watson, and F. Zachariasen, 
Sound Transmission Through a Fluctuating Ocean,
Cambridge University Press, Cambridge, 1979. 

FGPSbook
J.-P. Fouque, J. Garnier, G. Papanicolaou, and K. Sølna,
Wave Propagation and Time Reversal in Randomly Layered Media,
Springer, New York, 2007

garnier_evan
J. Garnier, 
The role of evanescent modes in randomly perturbed single-mode waveguides,
Discrete and Continuous Dynamical Systems B 8, 455–472  (2007). 

GP07
J. Garnier and G. Papanicolaou, 
Pulse propagation and time reversal in random waveguides,
SIAM J. Appl.  Math. 67,  1718–1739 (2007).

garrett
C. Garrett  and W. H. Munk, 
Internal waves in the ocean,
Annu. Rev. Fluid Mech. 11, 339–369 (1979).

gomez
C. Gomez,
Wave propagation in shallow-water acoustic random waveguides,
Commun. Math. Sci. 9,  81–125 (2011).

jensen
F. B. Jensen, W. A. Kuperman, M. B. Porter, and H. Schmidt, 
Computational Ocean Acoustics,
Springer, New York, 1993. Chap. 5.

kohler77
W. Kohler and G. Papanicolaou,
Wave propagation in  randomly inhomogeneous ocean,
in Lecture Notes in Physics, Vol. 70,
J. B. Keller and J. S. Papadakis, eds.,
Wave Propagation and Underwater Acoustics,
Springer Verlag, Berlin, 1977.






rb04
M. D. Richardson and K. B. Briggs,
Empirical predictions of seafloor properties based on remotely measured sediment impedance,
AIP Conference Proceedings 728, 12–21 (2004).


wilcox
C. Wilcox, Spectral analysis of the Pekeris operator in the theory of acoustic wave propagation
in shallow water, Arch. Rational Mech. Anal. 60, 259–300 (1976).








