Effect of the Acoustic Environment on Adjoint Sound Speed Inversions

The recent prevalence of low cost robotic platforms such as oceanographic gliders has increased the availability of long–term measurements of the ocean environment. Gliders can take direct measurements of the ocean sound speed environment, which is of interest in many ocean acoustic problems, including source localization and tomography. These measurements, however, have a low spatial–temporal resolution that makes them difficult to use directly. These measurements have the potential to provide an accurate environmental parameterization for acoustic inversions, which could in turn be used to measure the sound speed field at a much higher spatial–temporal resolution. This study uses glider measurements to provide the environmental parameterization used in the adjoint inversion method. The adjoint method calculates the gradient of a cost function describing the mismatch between observed data and acoustic model predictions with respect to the ocean sound speed. This gradient is a measure of how changing the sound speed at any point in the acoustic environment would affect this misfit. This cost function and its gradient information is then used as inputs to a numerical optimization routine, which efficiently finds a local minimum. There are two challenges of this method addressed in this study; the first is restricting the search space of this inversion. Proper parameterization of the inversion will ensure that the local minimum found in the numerical optimization routine is the correct result of the inversion. This parameterization allows for the combination of the relative strengths of both methods of measuring the sound speed field, the robust direct measurement of the glider and the near instantaneous result of an acoustic inversion. A covariance matrix is created from glider measurements of the range dependent sound speed field, which is then decomposed into an empirical orthogonal function (EOF) base. The mean profile and the significant EOF bases then form the search space of the adjoint method. The second issue is the proper treatment of the acoustic interaction between the ocean and its sea floor. The adjoint method uses the implicit finite difference form of the Parabolic Equation, which has a few possible bottom treatments. Two simple bottom interface treatments are the local boundary conditions of McDaniel and Lee [1] and the non–local boundary conditions of Papadakis et al. [2]. Local boundary conditions treat the interface by altering sound speed values at the interface to account for density interfaces. This interface treatment was selected over the more complete treatment of non–local boundary conditions, which treat the interface with a reflection coefficient. This decision was based on the relative simplicity of the testing environment, which did not require sophisticated bottom treatments. Finally, the performance of the adjoint method was tested by numerical simulation with a number of different test environments. The adjoint performance was most sensitive to the range of propagation, and relatively insensitive to other environmental parameters such as source depth and bottom depth. These results suggest that the adjoint inversion method will perform consistently in appropriate testing conditions.

olution that makes them difficult to use directly. These measurements have the potential to provide an accurate environmental parameterization for acoustic inversions, which could in turn be used to measure the sound speed field at a much higher spatial-temporal resolution.
This study uses glider measurements to provide the environmental parameterization used in the adjoint inversion method. The adjoint method calculates the gradient of a cost function describing the mismatch between observed data and acoustic model predictions with respect to the ocean sound speed. This gradient is a measure of how changing the sound speed at any point in the acoustic environment would affect this misfit. This cost function and its gradient information is then used as inputs to a numerical optimization routine, which efficiently finds a local minimum.
There are two challenges of this method addressed in this study; the first is restricting the search space of this inversion. Proper parameterization of the inversion will ensure that the local minimum found in the numerical optimization routine is the correct result of the inversion. This parameterization allows for the combination of the relative strengths of both methods of measuring the sound speed field, the robust direct measurement of the glider and the near instantaneous result of an acoustic inversion. A covariance matrix is created from glider measurements of the range dependent sound speed field, which is then decomposed into an empirical orthogonal function (EOF) base. The mean profile and the significant EOF bases then form the search space of the adjoint method.
The second issue is the proper treatment of the acoustic interaction between the ocean and its sea floor. The adjoint method uses the implicit finite difference form of the Parabolic Equation, which has a few possible bottom treatments. Two simple bottom interface treatments are the local boundary conditions of McDaniel and Lee [1] and the non-local boundary conditions of Papadakis et al. [2]. Local boundary conditions treat the interface by altering sound speed values at the interface to account for density interfaces. This interface treatment was selected over the more complete treatment of non-local boundary conditions, which treat the interface with a reflection coefficient. This decision was based on the relative simplicity of the testing environment, which did not require sophisticated bottom treatments.
Finally, the performance of the adjoint method was tested by numerical simulation with a number of different test environments. The adjoint performance was most sensitive to the range of propagation, and relatively insensitive to other environmental parameters such as source depth and bottom depth. These results suggest that the adjoint inversion method will perform consistently in appropriate testing conditions.   Inversion methods are designed to measure the properties of a physical system that are difficult or impossible to measure directly. Ocean sound speed, which is determined by the temperature, salinity and pressure of the seawater, is an example of such a system. In this case the number of measurements required to completely characterize the environment of an acoustic propagation is overwhelming for reasons of costs and logistics. Ocean acoustic inversion methods are therefore designed to measure the sound speed properties of the ocean by measuring an acoustic propagation that has passed through the environment. The adjoint method is such an inversion method, which seeks to determine which physical system best recreates an acoustic propagation measurement.
Knowledge of the sound speed of the ocean between an acoustic source and receiver is useful in several applications; two common examples are source localization and ocean acoustic tomography. Source localization is a problem that arises in passive SONAR systems, which attempt to determine the location of a ship or submarine by the noise radiated from the vessel. This problem is very sensitive to uncertainties in the acoustic environment, and solutions are more successful in well characterized environments [4].
Another example is ocean acoustic tomography, an important field in oceanography. These experiments are designed to measure the range-averaged temperature of the ocean. Acoustic inversions methods are used to solve for sound speed environment, which is then used to determine the temperature of the ocean. This technique requires a conversion from ocean sound speed to ocean temperature. The ocean parameters that determine sound speed are well characterized in literature [5] and sound speed is accepted as a sensitive indicator of ocean temperature [6]. The adjoint method is only one possible approach to the inversion of sound speed on small scales. A general discussion of inversion methods is required to consider the relative strengths of the adjoint method and its restrictions. Tarantola [7] breaks down the general inverse problems in three steps.
i Parameterization of the system: Discovery of a minimal set of model parameters whose values completely characterize the system (from a given point of view).
ii Forward modeling: discovery of the physical laws allowing us, for given values of the model parameters, to make predictions on the results of measurements on some observable parameters.
iii Inverse modeling: use of the actual results of some measurements of the observable parameters to infer the actual values of the model parameters.
The first step of the parameterization of the system is closely related to any previous knowledge of the system available to the inversion. The level of accuracy and relevance of this knowledge affect which approach is best suited for a given problem. In general, inversion problems attempt to integrate all of the available information about the system before solving for the system parameters. Therefore, when an independent ocean observing system is available to the inversion it is desirable to incorporate its measurements into the estimation of the system. Gliders are an example of a recently developed technology for ocean observation that fit into this networking approach. These autonomous underwater platforms use buoyancy changes and hydrodynamic shape to carry out zigzag motions between the surface and some depth into the ocean with a net horizontal displacement [8]. Nominal horizontal speed is about 0.5 m/s with spatial cycle periods depending on the programmed pitch and immersion depths [9]. Unlike profiling floats, glider motions are controllable to a degree determined by the strength of the current field. Their endurance allow long duration measurement of gross sound speed features local to a test area. However, they have limited speed in both the horizontal and vertical direction, which makes their measurements best suited for statistical, rather than direct, measurements of the sound speed field. This information is complimentary to an acoustic inversion survey, which can be thought of as taking much quicker samples of the sound speed field. Gliders instead provide a means to parameterize the physical system being measured with empirical data, increasing the success rate of the inversion process.
The second point in Tarantola's list is the forward model used in an inversion.
This organization is appropriate when performing an inversion, but less so when designing an inversion method. In the case of adjoint inversions the forward model is dictated by the choice of the inverse method. A general discussion of inverse methods is necessary to clarify the relationship between the forward model and the inverse method.
The inverse problem is ill conditioned; meaning the number of unknowns is larger than the number of measurements. This means that there is no single solution to the inversion problem. The strategy used to solve an inversion problem is therefore a design decision, made with the aim of finding the most physically likely solution. A common strategy is to create a least-squares cost function of the mismatch of the observed parameter and the prediction of a forward model. In this case the forward model can be chosen from a variety of ocean acoustic propagation models [4]. The goal of the inverse method is to minimize this cost function.  speed field has 5 possible parameters each with 20 possible values, a complete   search will require over 3 million function evaluations. Global search methods use   a semi-random search of the entire search space, always saving the best solutions, and will eventually cover the entire search domain. This approach is time consuming but has the advantage that it does not preferentially search one area over another. This means that global search methods provide robust searches in cost function that have spurious local minima, and one global minimum.
In contrast to global search methods, Hursky's adjoint method is a local optimization method. Specifically, this is a gradient-based local optimization.
Gradient-based optimizations can be thought of as proceeding down the slope of the cost function until all directions increase the cost. The specific strategy of the optimization is often more advanced than this steepest-decent algorithm, but similar in principal. While local optimization strategies can be applied to any function using numerical differencing, it is much more efficient when the gradient is analytically calculated [10]. The purpose of the adjoint method is simply to compute the gradient of the cost function [11].
Local optimization methods require many fewer cost function evaluations than The collection of appropriate sound speed data for the validation of this inverse method remains an outstanding challenge. Complete validation of the inversion method would require a sampling of the sound speed field at higher spatial and temporal resolutions than those provided by the acoustic inversion. This data would confirm the results could adequately characterize the sound speed field with the level of resolution used. Glider sound speed measurements are well suited for a preliminary study of the adjoint method however, because they provide long term sound speed information at low cost.
This work investigates a proposed method of characterizing a marine sound speed environment using a glider platform and an acoustic system. Glider data was collected during the REP11 experiment, performed in the Gulf of Taranto in September 2011. This data was employed to assess the performance of acoustic inversion when coarse glider observations along the acoustic path were available.

CHAPTER 2 Previous Work
The adjoint method for range dependent sound speed inversions was introduced to the field of ocean acoustics by Hursky et al. in 2004 [12]. In this in-

The adjoint method
The adjoint method was introduced to ocean acoustics by Hursky et al. in 2004 [12]. This method of computing the derivatives of cost functions had been used before this in the field of meteorology [11], geophysical inversion [13] and others. The adjoint method is part of a larger group of methods which compute Fréchet derivatives. The Fréchet derivative is the derivative of a function which maps one normed vector space to another, in this case sound speed to complex pressure amplitude [12]. The Fréchet derivative has been used extensively in ocean acoustic problems, where it is often called the Born-Fréchet kernel or the sensitivity kernel. Examples of it's use appear in travel time tomography [14], and angle of arrival tomography [15].

Implicit finite difference parabolic equation
The implicit finite difference parabolic equation (IFD-PE) was introduced by Lee et al. in 1981 [16]. This two-dimensional solution to the parabolic equation is based on the discretization of the propagation field into a grid of equally spaced depth and range steps. The spacing of the depth and range steps do not need to be equal to each other, however, and the range step is generally larger than the step in depth. In all cases the step size is a fraction of the wavelength of the acoustic signal, which for the case of the 400 Hz signal is approximately 3.8 meters.
This grid size is generally much smaller than the oceanographic features that are of interest in acoustic propagation. However, the computational cost due to this increase in model resolution compared with the earlier split step solution to the PE is offset in environments with significant range dependent features. This is because the split step approach requires the solution of a number of locally range independent sections, and becomes computationally inefficient which dealing with strong discontinuities in environmental parameters [4]. The IFD-PE handles these continuities gradually, with no increased computational cost.
Additionally, the angular restriction inherent to all parabolic equations can be relaxed in the IFD-PE using either the Claerbout or Greene square root approximation. The wide angle capability and good performance in range dependent environments made the IFD-PE a popular ocean acoustic propagation model in the 1980's, until it was eclipsed by the wide angle split step solution of Collins [17].

Interface treatment in the Parabolic Equation
When the IFD-PE was introduced by Lee et al. in 1981 [16], the treatment of the sound speed and density interface in bottom sediment layer was not addressed.
Since this paper, several different treatments of this interface have been introduced, each possessing relative merits and restrictions.
The first distinction between interface treatments of the IFD-PE is between local and non-local bottom boundary conditions (NLBC). The local bottom interface treatment, first introduced for horizontal interfaces by McDaniel and Lee [1], requires the numerical grid to extend into the bottom. Special treatment is given to the PE solver at the numerical mesh point on the bottom interface. This is a very simple treatment, but requires that the numerical grid containing the bottom to be at least the same size as the waveguide. Secondly, it requires that an artificial absorbing layer with high attenuation values be introduced at the bottom of this grid to prevent energy reflections caused by numerical effects of terminating this grid. Both of these effects require special attention when setting up the propagation environment, and also significantly increase computation cost.
Non-local boundary conditions, first introduced to the IFD by Papadakis et al. [2], introduces an impedance condition to the termination of the computational grid of the IFD. This impedance condition replaces the perfect pressure release boundary condition implied in the basic IFD formulation. The non-local impedance condition is computed to account for all of the sediment layers beneath the waveguide. This impedance condition is a more general treatment of the bottom effects than the local treatment, as it allows for the inclusion of shear effects in the sediment. Additionally, this treatment of the sediment interface was shown by Meyer and Hermand to have an analytic adjoint operator [18].
The fluid/fluid interface treatment required by the local boundary condition is an appropriate solution for low shear speed bottoms. This treatment was found by Tindle and Zhang to be valid for sheer speeds up to 200 m/s with small correction terms to the sound speed and density of the bottom [19]. This computational study will consider the effect of low sheer speed bottom, and so the local boundary conditions were considered sufficient to assess the performance of the adjoint method.

Range Dependent interface treatment
The treatment of range dependent bathymetric features in the PE has a long and involved history. The PE is attractive for range dependent features because they can be simply approximated by a stair step bottom. This approximation is simply implemented because the parabolic equation is a one-way equation. The depth of the interface can simply be changed at range indices that are designated as having a step.
This approximation was determined to produce erroneous results with bottoms with only mildly sloping bottoms when the Acoustic Society of America benchmark problem Wedge II was introduced by Jensen [3]. This was a deep to shallow wedge with an bottom which approximated sand, and had a rising slope angle of about 2.86 • . The stair step solution of the PE was found to agree with the one-way coupled mode solution, but diverged significantly from the full two-way solution.
The mechanism causing the difference in the two solutions was not backscattered sound, as expected. Instead, Porter et al. [20] demonstrated that the stair step treatment violated energy conservation with sloping bottoms.
The energy conservation issue arose from the vertical interfaces necessitated by the stair step approximation. While the horizontal interface treatment of McDaniel and Lee [1] and others could be shown to conserve energy, the vertical interface could not satisfy this condition in a one-way wave solution. The IFD-PE could be shown to match pressure at this interface, but not particle velocity. The split step PE solution could be shown to match reduced pressure, √ p at the interface. By testing a few possible matching schemes, Porter et al. [20] demonstrated that the simplest and most accurate matching scheme was to match p/ √ ρc at the interface.
Changes in the direction of development of the PE towards the split-step PE with Padé coefficients of Collins [17] mean the implementation of the p/ √ ρc matching was not widely discussed in PE literature. As a result, the simple pressure matching stair step implementation was used in this study. The implementation of an energy conserving stair step treatment of the irregular bottom is left for future work.

Parabolic Equation
The parabolic equation (PE) has found widespread use in the field of ocean acoustics over the past thirty years [4]. It's popularity stems from its ability to simply model propagation in range dependent environments. It is derived from a one-way wave approximation of the Helmholtz equation, and so it is has a few restrictions on it's use. Three major limitations are generally cited in literature.
These are: 1. the solution is limited to far field propagation, 2. that the range dependence of the problem is "weak", 3. in most solutions the angles of acoustic propagation have a limited aperture.
The first two limitation are covered in section 3.1.1, and the third is covered in section 3.1.2. Finally the IFD-PE solution of the PE is discussed in section 3.1.3.

Derivation
Linear acoustic propagation, which covers most problems of interest in ocean acoustics, is governed in general by the wave equation, a four dimensional equation.
The four independent variables are the three space dimensions, r = (x, y, z) and the time dimension, t. It is possible to reduce the dimension of this equation by modeling the acoustic source as a harmonic point source. This leads to the frequency domain wave equation, the Helmholtz equation, symbol parameter p(r, z) complex pressure field c(r, z) sound speed field γ(r, z) index of refraction squared (IORS) f frequency of acoustic signal p n range n pressure u n range n IORS N final range m experimental pressure H measurement operator Table 1. Common variables Table 1 is a list of definitions of the most commonly used variables.
In cartisian coordinates the operator ∇ is defined as In equation 1 the variable k is the wavenumber at the radial frequency ω, The frequency dependency of the Helmholtz equation can be removed by only solving for a single frequency at a time. The adjoint inversion performed in this study is only investigated at a single frequency. This approach is suggested in Hursky et al. [12], who found that wide band sources had little effect on the accuracy of the inversion result. For a single frequency, the ideal source is a pressure injection at a single point, modeled as a Dirac delta function, δ(r s − r).
The treatment of the acoustic source will be detailed later in this paper, and for simplicity only the homogeneous version of equation 1 will be considered. The inclusion of acoustic sources requires only minor modification of the derived result.
In many cases it is possible to reduce the dimensionality Helmholtz equation further by assuming that the acoustic propagation is only effected by the environment between the source and receivers. In many real ocean environments there are no significant cross path features which can refract sound from one plane of propagation to another, and this assumption is often used. In two-dimensional problems there are is a choice between two common coordinate systems, cylindrical and cartesian. In this discussion the cylindrical coordinate system is used, and all spatial variables are designated in (r, z). This choice of coordinates assumes a point acoustic source, and that the problem is radially symmetric.
Using the cylindrical coordinate representation of ∇ 2 , the homogeneous form of equation 1 is, To derive the one-way parabolic equation first express the acoustic pressure, φ(r, z) as two functions, The second function v(r) contains the range dependent spreading of acoustic pressure common to all ocean acoustic problems. When acoustic pressure is contained by the top and bottom boundaries acoustic pressure spreads cylindrical from the source, leading to an acoustic pressure reduction related to 1/ √ r. The function p(r, z) is the complex pressure field that contains pressure modifications specific to an acoustic environment.
Substitute this definition of φ(r, z) into equation 3 to get a separable differential equation, The reference wave number squared, k 2 0 , is used as the separation constant to solve equation 5. The left hand bracketed term leads to the expression of v(r) as v(r) = H (1) 0 (k 0 r), or a zeroth order Hankel function of the first kind. This function represents the range dependent pressure reduction due to cylindrical spreading [21].
The far field approximation of the Hankel function is defined in equation 6, and this formulation has a simply defined derivative.
This approximation allows for the simplification of the term (1/r + (2/v)v r ) to 2ik 0 . This approximation makes the parabolic equation unsuitable for near field solutions.
The The two way Helmholtz equation is reduced to a one way equation by factoring equation 7. Representing equation 7 using differential operators, Equation 9 may be factored into a forward and backward propagating wave components if the following identity holds, This is the communicator of the operators P and Q. This identity holds exactly in the case of a range independent sound speed field, n ≡ n(z). The parabolic equation assumes that the range dependence of the sound speed field is small enough that the non-zero communicator term is negligible. Using the identity of equation 10 and keeping only the outgoing wave term, equation 7 is factored as Equation 11 is an exact solution to the Helmholtz equation in the far field, for environments with no backscattering and range independent sound speed fields.
The backscattering and range independent requirements are violated frequently in practice, and relaxed to a requirement of "weak" range dependence [4].

Wide Angle Parabolic Equation
The final step of the derivation of the parabolic equation is to introduce an approximation of the differential operator, Q. This is commonly done by rewriting equation 11 to include the operator q, with Q = √ 1 + q. The one way parabolic equation is first rewritten to explicitly include the square root operator, This formulation allows for several different approximations of this operator.
Possible forms of these approximations are the Taylor series, a general rational function approximation, or a Padé series approximation. The choice of approximation limits the angle of acoustic propagation of the model.
The most basic approximation of Q is the first order Taylor series expansion suited for the adjoint method, and is investigated in this paper.
A higher order Taylor expansion requires the computation of higher order q terms, which can lead to computational issues [4]. Instead, other approximations based on the first order of q are used. The most basic of these is the rational function approximation.
This ration function approximation leads to the following wide angle parabolic equation, With the appropriate choice of coefficients this approximation is equivalent to the Trappert approximation, but other choices lead to more accurate approxima-

Implicit Finite Difference Parabolic Equation
The wide angle parabolic equation in 14 requires a definition of the inverse This can be done using a finite difference scheme. Lee and McDaniel [21] showed that an implicit finite difference scheme is necessary for the parabolic equation to be stable. The stability of a numerical solution states "the difference between the theoretical and numerical solutions remain bounded as the range step n increases, provided the range increment ∆r remain fixed for all space steps." [21]. The proof of the stability of the implicit finite difference scheme also showed that simpler explicit solutions of the parabolic equation based on a first order Taylor expansion were not stable.
The implicit finite difference scheme is based on a Crank-Nicolson scheme, The form of the PE shown in equation 15 is found using the central difference approximation of both the second derivative in depth and of the first derivative in range.
Factoring equation 15 to separate the pressure vectors by range step leads to the implicit equation This equation is then simplified algebraically.
To reach the final form of the IFD-PE requires the assumption that the operator q is constant over a range step. This is assumption means that the index of refraction can not vary with range with scales less than ∆r. This is a less restrictive assumption than the "weak" range dependence already required by the PE.
This requirement is reasonable because the size of the oceanographic features of interest are at least on the scale of the acoustic wavelength, λ, and the range step, ∆r, is a fraction of λ.
The complete form of the IFD-PE is quite long due to the number of points needed to construct a difference equation. This complete form is found in section 3.3.2, however, for much of our discussion a general matrix form of the IFD-PE is sufficient for understanding. This matrix form is shown in equation 17, The main diagonal of matrix B and C contain the index of refraction information. These matrices are also used to compute the second order depth derivative contained in the operator q. The matrices B and C are similar to each other, like the terms on both side of equation 16. The main diagonal is proportional to γ(r, z), the index of refraction squared (IORS). The conversion between sound speed and IORS is simply The existence of the inverse of matrix B means that the IFD-PE can be written in explicit form, with F n = B −1 n C n .
The parabolic equation is performed on a grid with constantly spaced range and depth vectors, r ∈ R r and z ∈ R z , respectively. The distance between grid points, ∆r and ∆z, is measured in meters, and is different for both the range and depth vectors. A vertical grid spacing of λ 4 and horizontal spacing of λ/2 was found to provide the limit of resolution of the result. The IORS field is then discretized to a grid γ ∈ R r×z for input in the parabolic equation.

Formulation of the Adjoint Method
The IFD-PE is simply defined as a vector matrix product which steps forward the acoustic pressure to the next range step with each multiplication. This

Tangent Linear Model
While the pressure field is only measured at one range, p R , changes in γ i,j at any point in the modeled area will effect this measurement. It will be necessary to have a linear model of this effect for all i and j before finding the cost function gradient using the adjoint method. Before addressing this non-linear problem, consider the linear case when the field γ is known, and the pressure field at all points in the PE grid is allowed to vary. The pressure at range p i would be defined as Using this definition with equation 19, results in The propagation of the perturbation in pressure at range i to range R is given in equation 22.
The second form of this equation introduces the notation F(R − 1 : i) = The linear perturbation analysis shown in equation 22 is exact if the matrix F is known for all ranges. This formulation does not match the reality of the experiment, however. First, there is only one source of pressure field in this experiment, and it is considered known. Second, the provision that F n is known for all n is equivalent to knowing γ, and thus eliminates the need to solve for the sound speed.
Perturbation analysis of equation 17 will show that perturbations in the IORS field lead to perturbations in the pressure field, and explain the sources of the pressure field perturbation seen in equation 21. This effect is non-linear, and requires a local tangent linear approximation. As a preliminary step, the range stepping form of the parabolic equation means it is convenient to work with perturbations in γ at each range step. This is accomplished by denoting the perturbation in γ at range step n as u n ∈ R z , where z is the number of depth bins used in the PE mesh and introducing a perturbation in u n of the first order in .
u n = u n + ũ n p n = p n + p n + . . .
Both B n and C n are diagonal matrices with the ith diagonal element equal to (k 2 0 /2)ũ n (z i ), though they have opposite sign.
Carry the multiplication in equation 24 through and rearrange the first order terms in a range stepping form, The matrix G is constructed from the complex pressure information for the zeroth order sound speed field. This means that each run of the forward model has two effects. This first is to calculate the pressure mismatch between the model and the measured data, which is used to determine a cost using equation 29. The second is to redefine the matrix G. This constant re-linearization of the gradient allows the optimization to handle small non-linearities in the cost function [12].
The perturbation analysis of equation 26 leads to a form that is very similar to the measurement sensitivity to a pressure field perturbation, shown in equation 22. This is most simply illustrated for the first range step. The starter field is assumed to be known exactly in this analysis, and so δp 0 = 0. With no pressure perturbation, δp 1 = −G 0 δu 0 . There is a similar relationship between every perturbation δu n with the perturbation of pressure at the next range step, δp n+1 . The error propagation presented in equation 22 can then be written in terms of δu i , These effects are cumulative because of the linearity of pressure perturbation propagation defined in equation 21, In this form, the sensitivity of p R to δγ is expressed as a sum of the sensitivity of p R to every δu. This solution is a linear solution for the sensitivity, and is applicable for small perturbations of δu.

Evaluation of cost function gradient
The least squares cost function is given in equation 29, Each value m i is a complex pressure field measurement taken at the ith hydrophone. The value p i is the pressure measurement predicted by the forward model for a given set of environmental parameters at the same index. The adjoint method will calculate the gradient of J with respect to γ. This method is based on properties of the inner product. The inner product is written x, y , which should be interpreted as x H y. The notation x H is used to denote the conjugate transpose of x, (x * ) (The notation x * is the complex conjugate of x and x is the transpose operator).
In the following discussion it is assumed that the hydrophones are ideal sensors, and the forward model is exact. In this case, the observation vector m ∈ R m is simply a sampling of the pressure vector p ∈ R n at m discrete points at range With these definitions equation 29 may be expressed as an inner product, The following properties of the inner product will be used in the adjoint method. First, the adjoint operator of A is defined as A H . The following identity can be shown to be a property of the inner product, Secondly, for small perturbations δu, the perturbation of J is related to the gradient of J with respect to u, ∇ u J, This is the form of the first variation, which is the core of the adjoint method.
If δJ can be manipulated so that δu is one vector of the inner product, the other vector must be ∇ u J. The order of the vectors is not important, since both vectors in equation 32 are real and a, x = x, a . From the definition of J, Substituting the definition of δp R from equation 28 leads to The adjoint operator introduced in equation 31 can be used to move matrices from one side of the inner product to the other. Additionally, the distributive property of inner products, x 1 + x 2 , y = x 1 , y + x 2 , y , can be used to remove the summation from inside of the inner product.
Using basic inner product algebra the IORS perturbation at each range, δu i , is With these definitions, equation 34 is concisely written This formulation clearly shows the large dimensionality of ∇ γ J, which has rz degrees of freedom.
With this definition, the cost function gradient of equation 36 is rewritten, The form of the gradient in equation 39 is amenable to numerical calculation.
First run the forward model to define λ R . Secondly, recursively define all values of λ until λ 1 . Once the value of λ i is defined, the gradient is defined for this range step i by multiplying G H i λ i . The transformation of measurements from the real diagonal track to a vertical profile is common in glider work. The first step in the glider data processing was to determine a common number of profiles taken on every transit. An interpolation grid of 5 profiles was found sufficient to encode the oceanographic variability along  Spatio-temporal variability is mostly concentrated in the surface layers and it can be described as a slight deepening of the thermocline and a reduction of the sound speed above it due to environmental warming.  Figure 2. Location of the glider measurement plotted with 100 meter bathymetric contour. Area shown in figure is the same as that of the blue box in figure 1

Bottom interface treatment
The sediment along the ocean bottom has a complicated effect on acoustic propagation, especially for low frequency propagation problems. The acoustic impedance, ρc, of the sediment in the ocean bottom is much closer to that of the ocean, which means that acoustic energy in one medium may penetrate into the other. This is in contrast to the air-ocean interface, which has a mismatch in acoustic impedance at least 4 orders of magnitude, virtually ensuring no energy transfer between interfaces.
Since acoustic energy propagates into the ocean bottom, the simplest way to model the bottom interaction is to extend the computational grid from the ocean into the sediment layer. An artificial attenuation layer is added at the bottom of this computational grid to ensure that no energy reflects from the artificial pressure release interface created by the termination of the computational grid. With this groundwork in place, the next step is to correctly model the acoustic propagation at the density and sound speed interface of the sea floor.
The simplest model of sea bottom is as a fluid, which does not support shear wave propagation. This approximation is applicable in many practical environments, and can be used for sediments that support shear sound propagation up to about 200 m/s using the equivalent fluid approximation of Tindle and Zhang [19].
In the case of a fluid-fluid interface, there are two conditions which are required at the interface: the continuity of pressure, equation 40 and the continuity of particle velocity, equation 41. for the values in medium 2. There are also two separate pressures in each medium, p 1 and p 2 , though they will be equated according to equation 40.
Following the derivation of the horizontal interface from Computational Ocean Acoustics [4], solve for the first depth derivative, p z by taking a Taylor series expansion of p n m−1 upon p n m . For brevity of notation, the pressure values without explicit grid indices are used to denote p n m . Instead, the notation p 1 is used to indicate the pressure in medium 1. Using these definitions, the expansion is, Solving for the second derivative, and substituting this result into 7 yields the derivative p z in the first medium, This result is used to satisfy equation 41. The derivation for p z in medium 2 yields an equivalent result. These results are then used to construct a partial differential equation for the acoustic pressure at an interface node. This is done by first equating p 1 = p 2 = p in equation 43 and it's equivalent in medium 2.
After this substitution, the two derivative expressions are multiplied by their respective ρ −1 and then equated. This leads to the representation of the far field Helmholtz equation 7 shown in equation 44, using the operator G to contain the finite difference terms.
The operator G is defined using the following definitions, The operator G is then simply, This is a general formulation for the IFD-PE, valid both at interfaces and at non-interface grid points when ρ 1 = ρ 2 = ρ and n 1 = n 2 = n. In this form the interface treatment is analogous to the two-way Helmholtz equation of 9. This equation is the factored into a one-way equation using the same assumption that the operators p r and G commute.
The difference between the finite difference formulation and the general formulation of the PE is that the rational-function operator approximation of equation interested reader is referred to either reference [4] or [21] for a complete derivation. and The terms of w are related to the rational function approximation of equation 13 as follows, where b is the index at the bottom. The bottom index of C is also the same.
Also, while the gradient of the cost function with respect to sound speed in also calculated for the bottom, it is simply set to 0 in this case, because this property is assumed to be known.

Dimensional reduction
Sometimes, phenomena which appear complex are actually be governed by a few simple variables. The fundamental assumption that justifies the dimension reduction is that the sample actually lies, at least approximately, on a space of smaller dimension than the data space. The goal of dimension reduction is to obtain a low-dimensional, compact representation of the data. This focuses efforts on the most representative portion of the dynamics, simplifying computations and, in general, increasing the accuracy of the results [22] by removing dynamics that are not observed from the search space.
Much of the expected uncertainty in the problem under consideration is associated with the lack of an accurate characterization of the environmental variability.
In this work, the variability of the sound speed is considered as resulting from a weakly stationary or second order Gaussian process defined by a covariance matrix. This covariance matrix can be estimated from the set of glider observations using the Maximum Likelihood Estimator (MLE).
The first step in defining the covariance matrix of the sound speed of the environment is to determine a measurable representation of this field. The measurements taken by a glider have a limited spatial sampling period, and a specific number of profile measurements are expected for a given experimental range. Each measurement of the sound speed field is then represented by a N sound speed profiles, u i , arranged in increasing range from the acoustic source. The full sound speed field used in the PE, γ is created using an interpolation routine, F .
The dimensionality of v is N z, or more simply v. An interpolation routine, F : R v → R rz is necessary because a measurement is not available for every value of γ. A simple linear interpolation is used in this study, but other interpolation methods may be substituted in F . The form of the function F is is not important to the minimization routine, which operates in the domain of the EOF bases, and not γ. For a given F , the knowledge of v is equivalent to knowing γ, and the vector v may also be referred to as the sound speed field.
The definition of the sound speed field in equation 53 is valid when considering either the direct sound speed measurement, or its perturbation,ũ. For clarity in the discussion of the covariance estimation, the complete sound speed profile defined in equation 23 is used. To calculate the covariance first define the sound speed detrended measurement matrix, M ∈ R v×n for n complete glider tracks, v i .
Using the spatial mean of the measurements,v ∈ R v to de-trend the measurement matrix.
The MLE of the covariance matrix of M is then given by  While equation 56 may be used to exactly recreate each measured profile, v i , these profiles may be estimated accurately using only a few principal EOF bases.
For discussion assume both the matrix D and V are sorted in decreasing magnitude of D. This behavior means that each succeeding EOF base is less significant than the the one before it in determining the final shape of the sound speed profile.
If the EOF bases are constructed from highly correlated data many vectors are of such small magnitudes that they may be considered numerical artifacts. The decomposition of the covariance matrix into its eigenvectors therefore provides a way to reduce the dimensionality of the initial problem by removing these numerical artifacts from equation 56 and projecting the dynamics onto the most representative EOFs.
It remains a problem, however, to determine which EOFs should be considered important. For a EOF to be considered significant in estimating the full shape of v it must achieve a reasonably large percentage of cumulative energy content.
This requirement can be made into a qualitative test by comparing the eigenvalues of the matrix D with those created by a white Gaussian noise (WGN) process.
A Bootstrap method with a null hypothesis given by WGN has been employed  After reducing the description of the sound speed field to a few principal components in the EOF approximation, the dimensionality of the gradient of the cost function may also be reduced. Introduce the chain rule of a gradient, The gradient ∇ α J ∈ R q relates changes in the magnitude of the EOF bases to changes in the cost function, J, equation 29. This is the final form of the gradient used in the adjoint inversion method. The matrix J α ∈ R nz×q is the Jacobian matrix, which relates a change in the magnitude the EOF bases to δc.
The Jacobian J c relates changes in γ to changes in sound speed. This second Jacobian is necessary because the adjoint method solves for the gradient of the cost function with respect to γ, where as the EOF decomposition is performed on sound speed.
The relationship between the IORS and sound speed was introduced in equation 18. The matrix J c is composed at the derivative of the sound speed at every index (r, z) with respect to the IORS at every index (r, z) as well. While this definition requires a nz × nz Jacobian matrix, there is a direct relationship between sound speed and IORS, and so the Jacobian only has nz non zero values.
The Jacobian J c is therefore a sparse diagonal matrix, and so J c

Error of EOF cutoff
The estimation of the total IORS field using EOF bases introduces a representation error component of the observation due to unresolved scales. This means that the uncertainty associated with an acoustic pressure measurement at sea may be larger than the error of the sensor. The present study is incapable of deterministically characterizing all the spatio-temporal scales of the environmental variability and it is therefore important to consider this representation error when finding an appropriate cost functions for well-conditioned inversion problems. The representation error expected on the sampled pressure field can be estimated given the dimensional reduction of the IORS field. This study illustrates the importance of this error term as an extreme case and assumes that there is no measurement noise at all. The variance caused by the EOF cutoff is the only source of noise and accounts for the maximum pressure resolution of the inversion model.
A Monte Carlo simulation of 100 trials is used to determine the maximum resolution of the inversion allowed by the EOF cutoff. First, create n = 100 profiles using full covariance matrix C, introduced in equation 55. Arrange these profiles in the matrix K ∈ R v×n , where each column in K is a simulated profile.
Second, project these profiles onto the restricted EOF space, E. The magnitudes of these projections are stored in the matrix A ∈ R n×q , where q was previously defined as the restricted number of EOF bases, The  Table 2. Measured pressure variance between IORS fields constructed using complete covariance statistics and their projection onto EOF bases.

Numerical Simulation of Adjoint Method
Numerical simulation is an useful method to determine the expected performance of the adjoint method in various situations.
The first step of a numerical simulation of an inversion is to simulate a possible sound speed field. This is done using the complete EOF basis measured by the gliders. The IFD-PE is then used with the simulated sound speed field to calculate the simulated measured pressure. This is the vector m used in the least squares cost function, equation 29.
The reference solution of EOF magnitudes is found by projecting the EOF bases onto the simulated sound speed field. These EOF magnitudes are denoted The result of the adjoint inversion can be directly compared with these reference solutions to determine its effectiveness. To measure the performance across many simulations, however, a single characterization of the inversion result is more useful. The root mean square error of the mismatch between the reference and in-

Irregular interface
The     The transmission loss for a reciever depth of 30 meters of the stair step IFD-PE is shown along with both the one-way and two-way COUPLE solution in figure 8.

Bottom interface treatment in the parabolic equation
As was the case with the range dependent propagation test, a mean offset between the two solutions was removed to account for different source normalizations between the models. While the one-way solution agrees well with the IFD-PE, there is a significant divergence with the two-way solution. The range dependent inter- face was not used for this reason, but should be valid when an energy conservation term is developed.

The effect of the enviornment on adjoint inversions
While the adjoint method can be expected to work in some acoustic environments, such as the one investigated by Hursky et al. [12], there need be limits on what conditions it works under. The limits of environmental conditions under which the adjoint method can be expected to converge to the correct result are explored in this section. The adjoint method was found to be highly sensitive to the range of the acoustic propagation.
When the adjoint method was introduced to problems of ocean acoustics by Hursky et al. [12], a very simple test setup was used to confirm that the method could correctly invert for the range dependent sound speed field. A 120-meter deep sound channel was bounded on both the top and bottom by a pressure release boundary condition. The transmission from a 400 Hz source was received at a distance 2 km away by a 32-element array that spanned the entire water column.
The basic test setup used in most computer simulations is shown in figure Modeling the bottom interface as a fluid/fluid interface with both a sound speed and density discontinuity increased the realism of the test setup. The first test with this new bottom interface treatment was to investigate its performance in a test case where the pressure release bottom treatment could correctly invert for the sound speed field. A 2 km inversion experiment was run with 5 EOF bases to compare the two treatments directly. The first 11 inversion runs used the sound speed fields measured by the gliders to simulate measured pressure. The rest of the inversions were performed on sound speed fields that had the same statistics of the measured fields, explained in section 3.3.3. The RMSE error, shown in equation 64, was used to compare the two results. This is shown in figure 10.
The performance of the adjoint method is comparable between the two bottom treatments. In general the pressure release bottom has slightly less RMSE than the results obtained using a sandy bottom, but both methods performed well in this simple test environment. There were no obviously divergent answers in this test scenario, with the largest error of a little more than 0.2. Compared to the mean RMSE value of about 0.1, this can be considered to be an acceptable answer.
The next set of experiments tested the range limit for which the adjoint method could be expected to converge to the correct result. A series of tests was run by expanding the range spacing between the five sound speed profiles shown in figure 9. The fields between these profiles were generated using an expanded linear interpolation routine.
There is a shift towards lower RMSE values in the 4 km results when compared with the 2 km results shown in figure 10. This is mostly because 9 EOF bases were used in this simulation study, instead of the 5 suggested by the null hypothesis test. This is discussed later in this section. Finally, a sandy bottom was used in both these tests. The result of range test is shown in figure 11. of scaled RMSE. The four cases where the RMSE is significantly larger than the mean RMSE value can be identified as outliers by inspection. This heuristic error criteria implies that 6 km is at the limit of the viability of the adjoint inversion method.
Finally, the RMSE criterion provides another insight into proper selection of the number of EOF bases to use in an inversion. The null hypothesis test showed that 4 EOF bases met a hard cutoff, and 5 could be used if this was stretched slightly. Increasing the number of EOF bases beyond this number continued to lower the RMSE value, however. Trial and error testing with increased numbers of EOF bases suggested that 9 EOF bases yielded a lower mean RMSE value over the suggested 5. Increasing the number of bases beyond 9 did not significantly change the results further. Finally, despite a mean decrease in RMSE value, the inversion results with more EOF bases still diverge on the same runs as the baseline case.
A comparison between inversion results with different numbers of EOFs is shown in figure 12.
The RMSE criterion for investigating inversion results is useful for observing the relative effects of environmental parameters. However, it cannot be viewed as an absolute test of the inversion performance. A more rigorous error analysis would provide a measure of the improvement to environmental understanding resulting from using the adjoint method. Identifying a specific application of the inversion results would provide a more objective rubric of performance. This would allow for the consideration of the penalty of divergent results, when the adjoint method converges to a minimum in the cost function that is unrelated to the minimum associated with the actual environmental parameters. In these cases the mean profile may be a better estimate of the environment than the adjoint result.
The performance of the adjoint method was shown to be largely unaffected by the introduction of bottom boundary conditions. This is an important step towards application of this method on experimental data. Secondly, simulation studies of the adjoint method using the RMSE error criterion show that the range of the environment strictly limits the adjoint method. Other possible parameter changes such as source depth, water depth and bottom type were investigated but did not show a significant effect on inversion performance. Within the proper operating conditions the result of the adjoint method is a consistent estimate of the actual sound speed environment, which is a useful measure of oceanographic properties and may improve the performance of other acoustic systems.

CHAPTER 5 Conclusions
The effects of a fluid/fluid interface were successfully included in the adjoint method. This addition to the method had minimal effects to the overall performance of the adjoint inversions. However, by including the bottom interface in the inverse method the adjoint method becomes applicable in more realistic propagation environments.
The adjoint method was shown to perform consistently in short-range environments with small range dependent variations in the sound speed field. It became unstable at longer ranges. In small propagation ranges it was seen to be relatively insensitive to the parameters of the test setup.
The experimental setup used to verify the adjoint method used a simulation study in which environmental statistics collected by gliders were used to simulate range dependent environments. This test setup was used to verify the basic viability of the adjoint method, and some of it's sensitivities to the environmental setup of the experiment.
However, the glider surveys should be considered very rough measurements of the sound speed field because of the limited spatial and temporal sampling frequency of these devices. A more desirable setup for testing the viability of realistic inversions based on glider surveys would use environmental statistics measured by much higher resolution techniques, such as a ship towed profiling device. The glider based adjoint inversions could be considered successful if the sound speed measurements taken by gliders were sufficient to resolve the more finely sampled data to the expected resolution.
The root-mean squared error criterion used in this study could show the relative performance of the adjoint inversions in different test setups. However, this error criterion does not provide a measure of the inversion performance that is relevant to any specific application. This heuristic criterion means that it is difficult to accurately say when the adjoint method becomes impractical or unproductive.
A more descriptive error measurement that could describe the cost of false results, perhaps designed for a specific application of the result, would provide a more rigorous error analysis.