Performance study of ray-based ocean acoustic tomography methods for estimating submesoscale variability in the upper ocean

Knowledge of variations in sound speed is essential for accurate predictions of sound propagation in the ocean, which is needed for various acoustic remote sensing applications.^1–3 Measurement and characterization of sound speed variability associated with submesoscale ocean features (with typical length scales of 100 m–10 km and time scales of hours to a few days^4–6) is currently an operational challenge, as such small-scale spatiotemporal variability cannot easily be accurately sampled, and often is aliased, by common oceanographic survey techniques (e.g., oceanographic floats). To address this challenge, we numerically investigated the estimation of range-dependent SSPs in the upper ocean over short ranges (< 5 km) using both the classical ray-based ocean acoustic tomography (OAT) method⁷ (with a fixed source-receivers configuration and fixed referenced sound speed) as well as iterative OAT and adaptive OAT formulations (i.e., when the sources and receivers configuration can evolve through time or space, which can be realized using mobile autonomous platforms) for comparison. To do so, a regional oceanographic simulation of the DeSoto Canyon circulation in the Gulf of Mexico is used to construct three-dimensional (3D) SSP variations spanning a winter month that exhibit significant submesoscale variability. For simplicity, hereafter three-dimensional effects such as horizontal refraction effects of the propagating acoustic rays are neglected and the simulated OAT inversions are conducted only in 2D vertical ocean slices.

In practice, sound speed variability occurs on a wide range of spatial and temporal scales and is driven by numerous types of ocean dynamical processes. For example, mesoscale flows in the form of eddies and fronts with typical dimensions of 10 to 200 km, time scales from months to years, and vertical velocities of 1–10 m/day, can spread and deflect ray paths.⁸ Furthermore, in the upper ocean, submesoscale features, again in the form of eddies and fronts, but with characteristics of length scales of 100 m–10 km, time scales of hours to a few days,⁴ and vertical speed reaching up to 100 m/day, are associated with intense density gradients that can further enhance sound speed variability. Unlike persistent mesoscale ocean fronts, submesoscale circulations modify the upper ocean structure on much shorter time scales and are driven by complex and momentary instabilities.^5,6 Therefore, submesoscale processes, which are especially significant in the tactical upper ocean (i.e., typically depths < 500 m), can affect underwater sound propagation in multiple ways, notably over (relatively) short propagation distances^9–12 due to the relatively small wavelengths of (most) man-made underwater sounds. For example, even at very low frequencies, e.g., 10 Hz, the corresponding acoustic wavelength in water is only $≈ 150$ m, which is on the smaller end of the conventional submesoscale range. Furthermore, ocean sound speed is typically estimated using empirical formulas based on local values of standard oceanographic parameters (such as temperature, salinity, and pressure)^1,2 and adequately measuring it as a function of the rapidly evolving submesoscale field would require a very dense spatial and temporal sampling of these oceanographic parameters throughout the volume of interest. Such a dense sampling is especially challenging using standard floats or moorings.

To alleviate this requirement, an alternative approach consists of directly leveraging the fundamental relationship between acoustic observations (e.g., arrival times measurements) and the actual volumetric spatiotemporal variability of the ocean sound speed profiles (SSPs). OAT, as originally introduced by Munk and Wunsch in 1979,⁷ refers to the process of using changes in acoustic measurements collected between multiple pairs of sources and receivers as inputs of an inverse problem aimed at reconstructing SSP variations with respect to a known reference environment within an ocean volume for 3D applications, or more commonly, for selected ocean slices for 2D applications. Specifically, the classical ray-based OAT method optimally combines the average variability of each ray path between each source-receiver pairs to estimate the sound speed anomalies (based on a reference sound speed model) within a selected ocean region. Much literature has been devoted to demonstrating the theoretical and experimental feasibility of OAT in various environments.^7,13–21 To date, however, most studies have focused on long-range (up to basin scales) OAT scenarios to primarily estimate mesoscale features, and very few studies have attempted to assess OAT performance to estimate the variability of sound speed driven by submesoscale processes.^18,20–23

One challenge associated with the “inverse problem”-based approach of OAT methods for estimating SSPs is that the relationship between acoustic observations (typically arrival-time measurements) and actual SSPs is fundamentally non-linear. However, fluctuations in fractional SSPs are typically small compared to baseline SSP values (usually $< 1 %$⁠). Hence, classical ray-based OAT methods simplify this approach by linearizing the relationship between the estimated SSPs perturbations and the measured variations in arrival times of stable rays propagating between the source and receiver pairs (the rays being identified by standard ray-tracing methods using the baseline reference environment, which is also referred to as the reference state). This provides a simpler and numerically tractable solution to the OAT inverse problem.^7,13–15 In addition, when solving the OAT problem, the reference state can be iteratively updated based on the estimated SSP anomalies of the previous (or initial) SSP values. This approach, known as iterative OAT^24–27 aims at further improving the estimated SSP variability with each iteration, but requires additional computations. However, it is a very tractable approach, given the numerical efficiency of the ray-based OAT method. Alternatively, the accuracy of the OAT solution can be further improved by constraining this inverse problem through additional regularization constraints of the sound speed values at specific locations. In practice, this strategy would correspond to using independent SSP measurements obtained from available CTD measurements at specific locations within the inversion domain, as described in Sec. III. Finally, it is also possible to gradually add additional source-receiver pairs at each iteration, which in turn increases the number of available ray travel-time measurements for the OAT inversion. This latter approach gradually enhances the effective ray coverage over the area of interest, which typically improves the accuracy of the estimated OAT solution and is referred to hereafter as adaptive OAT since the effective sampling of the inverse domain can be changed at each iteration. A practical strategy guiding this approach is to add acoustic or local SSP measurements in areas of the inversion domain where the OAT solution from the previous iteration appears to locally exhibit large SSP errors. In concept, these additional source-receiver pairs (and associated ray multipaths) in the upper ocean could be synthetically obtained using a suite of mobile autonomous platforms (or profiling)^18,20,28,29 equipped with a surface expression that allows for accurate positioning via a GPS link, as well as acoustic transponders and a CTD probe, as described in Sec. III.

To the best of our knowledge, an experimental data set that provides concurrent 3D SSP variations with sufficient temporal and spatial resolution to resolve submesoscale variability for testing OAT methods has not yet been available. Hereafter, we use a regional oceanographic circulation model covering the DeSoto Canyon in the Gulf of Mexico with a submesoscale permitting horizontal resolution of 0.5 km^30,31 and compute the corresponding 3D SSP variations calculated according to the TEOS-10 standard formula (using the Gibbs-SeaWater Oceanographic Toolbox)³² based on the standard modeled output variables (temperature, salinity, and pressure). We then compute the simulated ray arrival times for various source-receiver configurations distributed in a selected vertical slice of this ocean model using the bellhop ray tracing program.³³ The DeSoto Canyon area was selected as it includes significant density gradients resulting from the juxtaposition of relatively warm and salty waters advected from the equatorial Atlantic by the Loop Current, as well as the fresh riverine input by the Mississippi River System. Together, these density gradients contribute to the rich submesoscale field of fronts and eddies observed in all seasons.³⁴ Furthermore, the submesoscale circulation of this area has been well documented by drifter deployments following the Deepwater Horizon disaster (e.g., Refs. 35 and 36) and several modeling studies (e.g., including Refs. 34 and 37 and references therein).

In summary, the main objective of this study is to explore the feasibility of OAT methods to estimate the variability of sound speed in the upper ocean over short ranges (< 5 km) in the presence of energetic submesoscale ocean dynamics. A regional oceanographic simulation of the DeSoto Canyon circulation in the Gulf of Mexico is used to generate the 3D SSPs variations that are then used in the OAT simulations. The theoretical foundations for the various OAT formulations are presented in a consistent framework, allowing one to compare the estimated SSP values obtained through these OAT calculations for different source-receiver configurations and types of constraints. Therefore, the results of this study could be used to assist in the planning and design of future OAT experiments for comparable scenarios.

The remainder of this article is organized as follows. In Sec. II, the environmental parameters and source-receiver configurations used for numerical ray-tracing simulations and OAT inversion schemes are presented. Section III introduces the theoretical formalism and numerical implementation of the classical ray-based OAT method, as well as the iterative OAT and adaptive OAT variations. In Sec. IV, the OAT results are presented for the various selected configurations of environmental and acoustic sensors. Finally, a summary and discussion of the main results are given in Sec. V and closing remarks in Sec. VI conclude the study.

The DeSoto Canyon [see Fig. 1(a)] is located in the northern Gulf of Mexico and exhibits highly dynamical mesoscale and submesoscale variability throughout the year due to the proximity of relatively warm and salty waters advected by the Loop Current and associated eddies into the Gulf from the tropical Atlantic and the fresh and abundant riverine input by the Mississippi River System.^30,31,34 Notable SSP variability is found in the tactical upper ocean (i.e., depths < 500 m) where the submesoscale variability is greatest, especially nearby mesoscale eddies. This variability of SSP can affect underwater sound propagation in multiple ways,^9–12 and its estimation using OAT will be the focus of this study.

To this end, the ocean circulation in the DeSoto Canyon was modeled in the area shown in Fig. 1(a) (dashed box) using the Regional Ocean Modeling System (ROMS) in its Coastal and Regional Ocean Community model version (CROCO).³⁸ CROCO is a primitive-equation, free-surface, 3D terrain-following, split-explicit model and is configured, for this study, at the submesoscale permitting horizontal resolution of 0.5 km and 231 vertical levels with greater vertical resolution near the surface to capture the large density changes in the upper water column. The simulation is nested on a multi-year run at 1 km horizontal resolution covering the whole Gulf of Mexico, which was validated and described in further detail in previous studies (e.g., Ref. 39) The selected vertical resolution, while limited for a submesoscale permitting simulation, is comparable to what is commonly used by mesoscale permitting or resolving hindcast and forecast products, such as HYCOM.⁴⁰ The nested simulation covers 42 days starting on February 3, 2015. The first twelve days are discarded as a spin-up period, and the remaining winter days are used in this study. The modeled temperature, salinity and density fields are saved every 30 min and used to compute the 3D SSPs distributions using the TEOS-10 standard formula,³² creating T = 1439 time frames of sound speed environments.

To test the OAT methods, we selected a specific location in the southeast corner of the computational model, indicated by the red line in Figs. 1(a)–1(c). This location has been chosen partly because of the locally flat bathymetry that reduces potential terrain effects and partly because it is crossed by a diversity of mesoscale eddies and submesoscale features during the time window of interest. For simplicity, hereafter only the 2D sound propagation is modeled, and three-dimensional effects such as horizontal refraction effects of the propagating acoustic rays are neglected.¹ Therefore, the OAT inversions are performed only in 2D vertical ocean slices located at a constant latitude (27.5017°N) along the red line marked in Figs. 1(a) and 1(b). Specifically, ten 2D vertical slices of equal horizontal length (5 km), each separated by a distance of 0.5 km (equivalent to one range resolution cell of the ocean circulation model), were arbitrarily selected. They all had (nearly) the same reference SSP averaged over one month, as shown in Figs. 1(b) and 1(d), with a profile typical of a deep-water environment with an average SOFAR (sound fixing and ranging) channel axis centered at $≈$ 730 m. For all the ten vertical slices, most of the variability of the SSP, referred to as SSP anomalies in the context of the inverse problem of the OAT, is driven by submesoscale ocean characteristics and occurs in the upper ocean as expected [Fig. 1(c)]. In the remainder of this study, the first nine vertical slices provide only prior knowledge of the expected variability of SSP and are used to construct a database of SSP which, in turn, is used as a basis of empirical orthogonal functions (⁠ $EOFs$⁠) for dimensionality reduction and regularization of the OAT problem, as described in Sec. III. The last 2D vertical slice centered at a longitude of 86.3366° W is used to perform the OAT experiments. SSP discretization at depth, represented by the depth index z_i for the ith layer, is finer in the upper region with a grid size of 2 m from the surface to a depth of 100 m, followed by a grid size of 5 m from 105 to 1000 m, thus with a total of N_z = 231 points in depth for the discretized domain. The horizontal discretization of the SSP, represented by the index r_j is 500 m, where j = 1,…,N_r and N_r = 11. In the remainder of this study, the 2D discretized SSPs that cover this selected 2D slice are called sound speed maps (SSM). Given the selected vertical and horizontal discretization, the SSM to be estimated by OAT is composed of $J = N z × N r = 2541$ independent pixels, introducing the vectorized pixel index $j ∈ [ 1 , J ]$ representing all concatenated values of the 2D SSM.

The SSM $( z i , r j ; t k )$ in this dynamic range-dependent environment is thus a function of three variables: discretized depth z_i, discretized range r_j, and the discrete time step t_k used for ocean circulation simulations. For the sake of conciseness, the discretized subscripts will be removed later and the spatiotemporal dependence of the discretized SSM will simply be denoted as SSM $( z , r ; t )$⁠. To visualize the corresponding variability of SSM, a range-averaged anomaly δSSM $( z ; t )$ is computed for each time instant t in the 30-day simulation as the average of all N_r = 11 discrete ranges of the difference between the current SSM at time t and the monthly averaged reference SSM for this vertical slice. Figure 2(a) illustrates the temporal evolution of this range-averaged anomaly δSSM $( z ; t )$⁠: for example, for time 200 h $< t < 300$ h, the enhanced variability corresponds to an eddy that crosses the upper part of the water column (z < 200 m) during this time. To assess the performance of OAT in estimating SSP anomalies driven by submesoscale variability, visual inspection selected three different 7-daytime intervals during the month-long simulation to represent time windows with increasing overall magnitude of δSSM $( z ; t )$ anomalies. They are hereafter referred to as scenarios A, B, and C (see Fig. 2). Each of these three time windows was further divided into two parts: a 2-day training window (marked with subscript t) followed by a 5-day estimation window (marked with subscript e). This terminology is inspired by the development of machine learning methods for OAT.²¹ In this study, for each scenario A–C, the SSMs corresponding to the training window are only used to construct a $EOF$ basis (i.e., those training SSMs constitute prior knowledge of the expected variability of the SSM in the area). Furthermore, the actual SSMs from the estimation window are used to provide a ground truth for the estimated SSMs from OAT over this estimation window, as described in Sec. III. In general, these three scenarios exemplify different types of submesoscale features with various strengths and temporal evolution. It should be expected that the OAT method performs best when the variability of the SSM between the training window and the estimate window is the least. For example, scenario A is characterized by a training period that captures the spatiotemporal variability of the SSM which is comparable to the apparent variability of the SSM (primarily within [–1, 2] m/s) during the subsequent estimation period. However, the training period in scenario B appears to be slightly less representative of the SSM anomalies to estimate, which are within [–1, 3] m/s. Finally, scenario C shows the highest variability of SSM and thus the largest mismatch of SSM between the training and estimation periods, with higher anomalies in sound velocities within [–2.5, 4.5] m/s. This mismatch between training and estimation SSMs can be further quantified by computing the reconstruction error $δ c ̃$ root mean square error (RMSE) as shown in the next section and defined in more detail by Eq. (17).

The performance of OAT for the selected environment is investigated considering two main types of source-receiver configurations. The first setup shown in Fig. 3(a) is composed of a 20 source-element vertical line array (VLA) and a 20 receiver-element VLA with identical geometry, both spanning from 10 to 1000 m deep with an inter-element spacing of 52.1 m. This set-up represents a dense acoustic sensing geometry which provides a relatively dense ray coverage of the SSM pixels [see Fig. 10(a)] and is thus expected to yield the most accurate OAT estimates. The simulated ray arrival times for this selected source-receiver configuration are then computed using the Bellhop ray tracing program.³³ For each source-receiver pair, only the first two arrival times (direct and surface bounce paths) are used for the OAT inversion hereafter [Fig. 2(b)]. Other multipath arrivals (e.g., bottom bounces) are purposefully ignored, as in practice those later arrivals carry additional uncertainty (e.g., spans the water column from 245 to 695 m and its elements are spaced 50 m apart due to the typical uncertainty of most bathymetry databases) which exceeds the accuracy requirements for OAT. Hence, for this first setup, a maximum number of I = 800 arrival time measurements is possible. The second setup used in this study is a non-ideal—for OAT purposes—sparse distributed geometry [Fig. 3(b)] composed of three shallow sources positioned at 114, 322, and 531 m and a 10 receiver element VLA positioned from 62 to 1000 m with a constant 104 m element spacing distance, yielding a potential number of I = 60 arrival time measurements. This second configuration represents a more realistic and sparse acoustic setup, which could potentially be implemented in a synthetic fashion using few profiling sources and receivers. The adaptive methods presented in Sec. III use this second setup as an initial state, but then consider the successive addition of extra VLAs composed of identical transponder geometry (i.e., where each element can act as both a source or receiver in a round-robin fashion¹⁷) to sequentially augment the total number of arrival-time measurements for OAT acoustic. Each additional VLA provides 10 new sources and receivers spanning depths from 245 to 695 m with a 50 m element spacing distance. Hereafter we will consider up to three additional VLAs placed, respectively, at ranges r = 1500, 2500, and 3500 m. These extra VLAs could be implemented in practice by using a profiling transponder and relying on a “frozen ocean” assumption during the time required to collect the acoustics measurements between all the selected source-receiver pairs by profiling accordingly along the specified VLA aperture.

The primary focus of this study is to evaluate the performance of ray-based OAT methods to estimate submesoscale SSM variability using some acoustic sensing geometries shown in Fig. 3. However, the OAT estimates shown in Sec. IV represent only lower-bound errors for the various types of submesoscale variability present in the selected environment. In this initial performance study, a wide variety of error sources are not modeled, but they would need to be taken into account to fully quantify the uncertainty of OAT estimates, as extensively discussed in the literature.^7,14 Such errors would notably be caused, among others, by transducer positioning errors and unaccounted motion (e.g., array tilt due to tides or currents), imperfect time synchronization between the various source-receiver pairs, the impact of noise and available acoustic bandwidth on the accuracy of the measured ray travel times, and variability at scales smaller than those simulated by CROCO.

In general, the accuracy of the estimated SSMs from OAT depends on the selected parameterization of the simulated SSMs and, in particular, the number of time samples N_t (e.g., the number of days) used to form the aforementioned training set (which forms the $EOF$ basis), as well as the number $n Φ$ of EOF coefficients, which are ultimately kept to perform the expansion given by Eq. (6). Various empirical metrics can be used to select the optimal parameters N_t and $n Φ$⁠. For this study, these parameters were selected to minimize, in the test window, the error of the OAT estimates obtained using the methodology presented in Secs. III B and III C. Using this approach, it was found that an adequate basis $EOF$ could be calculated using only a 2-day training window (that is, N_t = 48 samples for each of the 10 slices shown in Fig. 1) for each scenario A–C (see Fig. 2). Furthermore, the following values for the parameter $n Φ$ for scenarios A, B, and C of, respectively, 27, 10, and 17 were found to be adequate to represent the variability of the sound speed for the selected environment, as well as to minimize, on average, the errors of the OAT estimations. However, it should be noted that the selected values for the parameters N_t and $n Φ$ generally depend on the environmental variability of interest. Furthermore, Fig. 4 illustrates in a simple way how changing either the parameter N_t or the parameter $n Φ$ can influence the parameterization of the sound speed variability for this study. For example, while keeping the duration of the training window at two days, Fig. 4(a) shows the variations of the magnitude of the singular values of the EOF basis matrix $Φ$ for the three scenarios A–C. The rapid decay of the magnitude of singular values confirms that most of the variability in SSM is well captured using less than 30 EOF coefficients for the three scenarios A–C, and thus the selected value of the parameters $n Φ$ mentioned above is consistent with the variability exhibited by scenarios A–C. Further discussion of Fig. 4 is given in the results in Sec. IV after having defined the OAT estimation methodology and the associated error metrics in Secs. III B and III C. Finally, it can be noted that alternate parameterizations of sound speed fluctuations have been proposed in the literature, notably using a set of horizontal and vertical fixed basis functions.^10,20,26 A potential benefit to this later approach may be the ability to parameterize a larger set of SSMs, which may otherwise not be accurately described using Eq. (6). An in-depth analysis of the effect of the selected parameterization of the sound speed fluctuations on the performance of the OAT methods when estimating submesoscale variability in the upper ocean requires further investigation in future studies.

Finally, we practically implement the solutions given by Eqs. (8) or (12), the optimal values for the penalty parameter λ (or λ₁, λ₂ for the multi-objective case) as well as the number of coefficients $n Φ$ for the EOF basis must be selected. In this study, for each of the three tested scenarios A–C, these parameters were selected as the optimal ones that minimize the RMSE of the estimation error made by the OAT, as defined in the Eq. (15), over the whole estimation window. Those optimal parameters were found to also minimize the error in residual ray travel times computed from ray tracing using either the a priori known SSM or the estimated SSM, thus confirming that the selected inversion parameters could also be determined in this alternative manner.^43,44

The true sound-speed anomaly $δ c t , j = c t , j − c 0 , j$ is the quantity of interest to be estimated in this study [see Eq. (1)], where $c 0 , j$ represent the reference (average) vectorized SSM and $c t , j$ represented the true vectorized SSM.

When written in terms of actual range-depth pixels of the SSM, instead of vectorized pixels, both the true sound speed anomaly $δ c t , z , r$ and the estimated sound speed anomaly $δ c ̂ t , z , r$ are each functions of three variables: discretized depth z, discretized range r, and the discrete time t used for the ocean circulation simulations. Hence, additional granularity for quantifying the spatial variability of the RMSE associated with which each OAT estimates can be obtained, for instance, by performing a spatial averaging across depth or range only instead of overall pixels as shown in Eqs. (14) or (15) or averaging over time only. For the sake of completeness, all additional error variables used in Sec. IV are individually defined below. The general notation convention used in this paper is that for each defined RMS error, the subscript indices $( t , r$ or $z )$ denote which of the physical variables is kept as a free variable while all other physical variables are being averaged over their respective total number of samples (i.e., T  $= 240$⁠, N_r = 11, or N_z = 141).

Next, we discuss the implementation of iterative OAT formulations, for which the reference state for the OAT inversion is iteratively updated based on the previously estimated solution, as well as adaptive OAT formulations, in which the sources and receivers configuration can evolve through time or space in between each successive iteration, for the environmental and sensing configurations previously discussed (Figs. 1–3).

The iterative OAT method²⁴ aims at improving the accuracy of SSM estimates, especially when the reference state (that is, the initial vectorized SSM map $c 0$⁠) significantly differs from the true SSM to be estimated. The procedure is inspired by the Levenberg-Marquardt algorithm,⁴⁵ which approximates the solution of a least-square problem by iteratively linearizing the forward model based on the previous solution found. Several studies have investigated the iterative OAT method in the context of long-range scenarios,^24–27 and have shown that the iterative OAT scheme generally improves the accuracy of the SSM. Here, we apply the same procedure to assess the accuracy of the estimated SSM in the presence of submesoscale variability over short propagation ranges. Specifically, the iterative OAT method first initializes the inverse problem variables with the reference state ${ τ 0 , E 0 , c 0 }$ associated with the iteration index n = 0, and then produces an estimated SSM $c ̂ n = 1 = c ̂$ based on one inversion using any of the OAT methods described in Sec. III. For each subsequent iteration $n ≥ 1$⁠, the algorithm updates this reference state (see Fig. 5) by using the estimated SSM $c ̂ n$ obtained from the previous iteration as a reference SSM and computes an updated set of ray arrival times $τ n$ between the selected source-receivers pairs as well as the corresponding linearized forward model $E n$ [using Eq. (5)]. Additionally, at each iteration both the optimal values for the penalty parameter λ (or λ₁, λ₂ for the multi-objective case) as defined by Eq. (8) or (12), and the number of coefficients $n Φ$ selected for the EOF basis can be updated.

The algorithm for the adaptive OAT method is a generalization of the one for the iterative OAT method, in which the acoustic sensing setup (that is, the specific location and number of selected source-receiver pairs) is allowed to change in between successive iterations. In this study, the adaptive OAT method is implemented using the sensor geometry shown in Fig. 3(b): for the initial estimation (n = 0) the acoustic sensing setup only uses a sparse combination of three shallow sources (positioned at 114, 322, and 531 m) and a 10 receiver-element VLA (positioned from 62 to 1000 m with a constant spacing distance of 104 m of elements). For each subsequent iteration $n ≥ 1$⁠, an additional 10 element VLA, where each element can act as both a source or receiver in a round-robin fashion,¹⁷ increases the acoustic sensing geometry and therefore improves the effective ray coverage for the OAT inversion. In this study we consider up to three iterations and the VLA added after each iteration has the same geometry but is placed at a different range r = 1500, 2500, 3500 m. Figure 3(b) displays the final sensing geometry for the third (n = 3) OAT iteration which now includes all three additional VLAs.

In the most general case, these added VLAs for the adaptive method can be added at an arbitrary location across each iteration. In practice, those added VLAs, or any desired acoustic transponder configuration, could be placed in regions of the domain that exhibit high SSM anomalies or variability to help better constrain the OAT inversion. However, for this study, we arbitrarily chose to keep the added VLAs at the same locations to facilitate the quantitative comparison of the OAT estimation performance between either the 3 × 10 or 5 × 20 initial sensor configuration, as shown in Secs. IV B and IV C to keep the overall sensing configuration as consistent as possible.

Figure 6 compares the accuracy of the various linear estimators given in Sec. III B when applied to scenarios A–C [see Fig. 2(a)] using the dense acoustic sensing setup of 20 sources and receivers [see Fig. 3(a)]. As explained at the end of Sec. III A, the estimators' hyperparameters $n Φ$ and λ (equivalent to λ₁ and λ₂) are chosen such that the inverse RMSE for the fully constrained adaptive method is minimized over the estimation window. Therefore, the following values of the parameter $n Φ$ for scenarios A, B, and C were selected, respectively, as 27, 10, and 17. Similarly, the values of the selected values for the hyperparameters λ₁ and λ₂, for scenarios A, B, and C were selected as [ $10 − 9$⁠; $10 − 8$], [ $10 − 9$⁠; $10 − 13$], and [ $10 − 10$⁠; $10 − 9$]. Figures in the first column [Figs. 6(a), 6(e), and 6(i)] show the temporal evolution of the range-averaged SSM anomalies $δ c t , p$ [see Eq. (22)] for each scenario. The figures in the second column [Figs. 6(b), 6(f), and 6(j)] show the equivalent estimations $δ c ̂ t , p$ obtained using the constrained Tikhonov estimator [see Eq. (10)]. The third column figures [Figs. 6(c), 6(g), and 6(k)] present the depth-dependent root mean square anomaly [see Eq. (23)] (green lines) and the estimation $δ c ̂ p$ RMSE [see Eq. (25)] for the different estimators (tik, cov, const. tik, and const. cov) for each method (other colored lines). Finally, the fourth column figures [Figs. 6(d), 6(h), and 6(l)] detail the anomaly time evolution RMS [see Eq. (14)] (green lines) and the associated estimation $δ c ̂ t$ RMSE for all estimators [see Eq. (15)] (other colored lines).

The temporal evolution of the reconstruction error $δ c ̃ t$ RMSE [as defined in Eq. (17)] during the test periods for scenarios A, B, and C are shown by the black dotted line in Figs. 6(d), 6(h), and 6(l). Overall, the average values of $δ c ̃ t$ RMSE across all estimation periods (⁠ $t = 1 : T$⁠) for scenarios A, B, and C were found to be 0.01, 0.14, and 0.41 m/s, respectively. Those RMSE values confirm that scenarios A, B, and C are expected to have an increased level of SSM mismatch between the training and estimation periods, which leads to a higher estimation error for the OAT inverse problem. As a point of comparison, Fig. 4(b) shows the temporal evolution of the reconstruction error $δ c ̃ t$ RMSE [see Eq. (17)] for the most variable scenario C for increasing duration of the training set from 1 to 15 days (while keeping the parameter of the number of EOF coefficients parameter $n Φ = 17$⁠). The displayed reconstruction error values indicate that changing the duration of the training period has little effect on the parameterization of the sound speed variability, even for the highly variable scenario C, for the first two-thirds of the test period. As expected, a gradual increase in the reconstruction error still occurs, and it can be noted that for the last third of the test period, using a longer training period (or more than 10 days) could slightly lower the reconstruction error. However, for convenience, the duration of the training period is kept constant to two days for all three scenarios A–C. Furthermore, Fig. 4 shows, as an illustration, the spatial variations of the reconstruction error [see Eq. (16)] for both scenarios A (c) and C (d) at the end of the test period (which represents the worst-case scenario), as the EOF basis obtained from the prior 2-day training period is least likely to be accurate at the end of the 5-day test period.

Based on the hyperparameters listed above, the averaged value of the OAT estimation errors increases as expected from scenarios A (⁠ $≈ 0.3$ m/s), to B (< 0.5 m/s) and C (> 0.5 m/s) proportionally to the increase in amplitude of the actual SSM anomalies across these three scenarios. The results for each scenario and the selected linear OAT estimator are summarized in Table I. It can be seen that for both the unconstrained [Eq. (8)] and constrained OAT estimators [Eq. (12)], using either the Tikhonov-based regularization or covariance-based regularization yields a nearly identical performance. Hence, in the remainder of this article, only the Tikhonov-based regularization is adopted. Results from this section confirm that most of the submesoscale features occurring in the 1 km upper-ocean channel have been correctly estimated in the ideal acoustic setup of 20 sources and 20 receivers [Fig. 3(a)], which offers a sufficiently dense coverage for the 5 km separation range, when applying the $EOF$ decomposition of the training set as spatial filter on the sound speed anomalies to be estimated. This corroborates the suitability of the OAT methods for this environment, provided that a sufficiently dense ray coverage and prior knowledge of the expected SSM fluctuations are available.

The results for the iterative and adaptive OAT methods (see Fig. 5) are presented for scenarios A and C for two different initial sensor configurations with different sensor densities. The first one is a 3 × 10 sensing configuration [as shown in Fig. 3(b)] which initially offers the lowest ray-path coverage of the inversion domain. The second one is a 5 × 20 configuration, where all 20 available receivers of the dense setup studied in Sec. IV A [as shown on Fig. 3(a)] are kept but only 5 out of the 20 sources are kept, still spanning the whole depth of the inversion domain down to 1000 m. These 5 sources are at specific depths: 10, 218, 479, 740, and 1000 m using. This second configuration of 5 × 20 initially offers an improved ray path coverage of the inversion domain compared to the configuration of 3 × 10.

Each iterative or adaptive OAT method is iterated four times (⁠ $n = [ 0 , 1 , 2 , 3 ]$⁠, including the initial estimation based on the reference state $c 0$ for the SSM). For the adaptive OAT method, the VLAs added at each iteration are identical and located as shown in Fig. 3(b) for either the 3 × 10 or 5 × 20 initial sensor configuration. For the iterative OAT method, no additional VLAs are added and only the 5 × 20 configuration or the 3 × 10 configuration is used for all iterations.

The green curves in Figs. 7(a) and 7(f) show the spatially averaged and time-dependent values of the initial RMS SSMs $δ c t$ [see Eq. (14)]. On the same figures, the yellow curves represent the initial estimation errors $δ c ̂ t | n = 0$ RMSE with respect to the tik. estimator [see Eq. (9)] for the 3 × 10 sensor geometry [see Fig. 3(b)]. The dashed lines represent the final estimation errors after the last iteration $δ c ̂ t | n = 3$ RMSE, for the iterative method (blue line) and the adaptive method (red line). Figures 7(b) and 7(g) show the ground-truth SSM anomaly $δ c t$ at the end of the 5-day estimation period (t = 696) for both scenarios A and C. Figures 7(c)–7(e) and Figs. 7(h)–7(j) show the inverse mapping error maps $δ c ̂ t | n$⁠, for the same time t = 696 and $n ∈ [ 1 , 3 ]$ [Eq. (26)] for the adaptive OAT method using const. tik. estimator. Note that the final inverse mapping error maps [obtained with Eq. (26) where n = 3] correspond to Figs. 7(e) and 7(j). Overall, these 2D error maps show a gradual improvement in the accuracy of the OAT estimates after each iteration for the selected configuration.

Interestingly, using the unconstrained OAT estimator [Eq. (8)] for the iterative OAT method (blue dashed lines) yielded slightly better performance than using the constrained OAT estimators [Eq. (12)] (not shown here for the sake of conciseness). On the other hand, the reverse trend was found for the adaptive OAT method (blue dashed lines): that is, the constrained OAT estimators [Eq. (12)], performed best. It is hypothesized that this may be due to a sensing configuration effect, as—by choice for practical considerations—the sound velocity constraints are only imposed where acoustic transponder elements exist using the formulation given by Eq. (12). For the case of the iterative OAT method, since the sensing configuration is fixed across iterations and remains sparse and equal to the initial 3 × 10 geometry, adding constraints only at a few discrete locations on the edges of the inversion domain leads to more estimation error. In contrast, for the adaptive OAT case, the sensing configuration evolves over iterations and spreads over the entire inversion domain as the number of iterations increases [Fig. 3(b)]. Hence, it gradually yields more distributed sets of constrained sound velocities locations across the inversion domain, which improves the accuracy of the estimated SSMs for an increasing number of iterations. Figure 8 further quantifies the influence on the types of constraints provided by the additional sensors of the three added VLAs (introduced after each iteration) on the precision of adaptive OAT estimators for scenario A or C (scenario B is omitted for simplicity). To do so, three different cases were simulated. The first case (blue line) corresponds to added sensors only providing local sound speed values (and no additional acoustic travel-time information); this would correspond in practice to the case of these sensors being only CTD probes and not actual acoustic transponders. The second case (red line) is the reverse, that is, the added sensors only provide additional acoustic travel-time information but no local sound speed values; this would correspond in practice to the case where these added sensors are only acoustic transponders [i.e., by setting the parameter $λ 2 = 0$ in Eq. (12)]. Finally, the third case (purple line) is the same as the fully constrained adaptive OAT case shown in Fig. 7, that is, when the added sensors provide both additional information on the acoustic travel time and the local sound speed values. As expected, Fig. 8 for scenarios A and C indicates the precision of the constrained adaptive OAT estimator [i.e., using Eq. (12)] is best when both constraints of acoustic travel times and sound speed values are used to perform OAT, especially for the highly variable scenario C. Thus, in summary, for the remainder of this study, all iterative results of the OAT method are obtained with the non-constrained Tikhonov estimator, while all adaptive results of the OAT method are obtained with the constrained Tikhonov estimator (for both the initial sparse 3 × 10 sensing configuration and the intermediate 5 × 20 sensing configuration).

On average, the relative improvement gain obtained by the iterative and adaptive methods is close to $20 %$⁠. The iterative OAT method improves more than the adaptive OAT method for scenario A. On the other hand, the two methods are almost equivalent in scenario B and the adaptive OAT method is the most effective with respect to scenario C. Overall, these results indicate a dependence between environmental variability, ray coverage, and the best estimator choice, as discussed in Sec. V.

A simple comparison is conducted to investigate to what extent the results for the classical OAT method using the dense 20 × 20 sensor configuration for the acoustic setup [as shown in Fig. 3(a)] can be outperformed (or not) by applying the adaptive OAT method using an initial 3 × 10 sparse configuration of the acoustic setup [shown in Fig. 3(b)]. Here, the const. tik. OAT estimator [Eq. (12)] is used for all cases. No additional iteration is performed for the dense 20 × 20 sensor configuration as the estimation errors are already very low after the first (i.e., n = 0) inversion (see Fig. 6). The results of this configuration using the classical OAT estimation are adopted as a baseline for comparing the performance of adaptive OAT estimations.

Figure 10 compares the temporal evolution of the RMSE $δ c ̂ t$ [Eq. (15)] of the estimated SSM obtained using either the dense 20 × 20 acoustic setup (blue line) or the initial (i.e., n = 0) sparse 3 × 10 acoustic setup (red line) or the final configuration (yellow line) for the adaptive setup [i.e., when all 3 extra VLAs have been added to the setup as shown in Fig. 3(b)]. For scenario A, a single inversion using the dense 20 × 20 setup yields the lowest estimation error. However, the trend is reversed for scenario B, and even more notably for scenario C, where the adaptive OAT estimations lead to the lowest estimation error. These results indicate that the error reduction obtained from the adaptive OAT method is more noticeable when the initial SSM anomalies [see Eq. (26)] are large, as for scenario C (see Fig. 2).

The results for the three cases shown in Fig. 10 can be further explained by considering the ray coverage maps $R cov$⁠, which show the number of eigenrays having crossed each pixel of the SSM [Fig. 11(a)–11(c)], since the effective ray coverage inherently drives the performance of the OAT inversion. Essentially, the classical dense setup (20 × 20) has fairly good ray coverage (i.e., over 100 eigenrays per pixel) over most of the inversion domain [Fig. 11(a)]. Note that because of refraction, the bottom center part of the inversion domain shows less dense ray coverage. On the other hand, the ray coverage for the distributed setup (3 × 10 initial) is initially very low (<50 eigenray per pixel) for most of the inversion domain [Fig. 11(b)], but is significantly enhanced (>100 eigenray per pixel over most of the domain) in the final configuration of the adaptive setup [i.e., when all 3 extra VLAs have been added; see Fig. 3(b)]. This improvement in ray coverage in the adaptive case explains most of the reduction in the time-averaged (over the 5-day estimation period) inverse mapping error $δ c ̂ ¯ | n$ [see Eq. (27)] between the values obtained for the initial estimation [n = 0, Figs. 10(e) and 10(h)] and the values obtained for the final iteration [n = 3, Figs. 10(f) and 10(i)], which are now comparable to the values obtained from the dense (20 × 20) acoustic set-up [n = 0, Figs. 10(d) and 10(g)] especially in scenario C (last row).

This study evaluated the performance of various OAT methods in the presence of submesoscale circulations. A summary and discussion of the main results are given below.

First, the results of this study show that since ocean submesoscale variability occurs on small spatial (100 m to a few km) and temporal (hours to days) scales, accurately estimating this variability with OAT requires both fairly dense acoustic sampling of the domain of interest (which is consistent with findings from previous OAT studies¹⁴) as well as selection of the reference point used to create the linearized forward model as close as possible to the onset of the estimation time period (which is key to provide accurate OAT estimations). Nevertheless, OAT can still be a viable alternative to conventional oceanographic sampling methods (e.g., local point measurements on moored or drifting platforms) given the high density of measurements required to capture this submesoscale variability over small spatial scales. Additionally, the EOF basis used to regularize the inversion problem should be selected to best represent the most current expected variability in the operational area of interest (e.g., a 2-day interval was used to construct the EOF basis in this study).

Second, the performance of the OAT estimation depends on the spatial extent and magnitude of the sound speed anomalies (see Fig. 6). As expected, large (in dimensions) and weak SSM anomalies are easier to reconstruct (scenario A). Intermediately, large features with relatively stronger SSM anomalies can be reconstructed with as much precision as small-sized anomalies with weak magnitude (scenario B). Last, highly variable, small, and strong SMM anomalies prove to be the most challenging to estimate (scenario C). In particular, the magnitude of the reconstruction error (that is, how accurately the SSM to be estimated can be described by the selected EOF basis for regularization purposes) is highly linked to the estimation error obtained from OAT. For the configurations tested in this study, a RMSE difference between the reconstruction error and the estimation error below 0.1 m/s appears sufficient to consider the estimated SSMs as near optimal.

Third, iterative and adaptive OAT methods offer the potential to improve the accuracy of the estimated SSMs for the studied environment in comparison to classical OAT methods (i.e., methods which consist of only a single inversion). Fundamentally, a sufficiently dense ray coverage of the ocean domain of interest is key to accurately estimate the rapidly evolving submesoscale features using any OAT estimator. Results from Figs. 10 and 11 show that the classical OAT method is sufficient if the magnitudes of the sound speed anomalies are weak (< 1 m/s) to moderate (⁠ $≈$ 1 m/s) (e.g., scenario A or B) and the ray coverage is fairly dense [e.g., see Fig. 3(a)]. The iterative OAT method, on the other hand, is most beneficial in the presence of low-to-medium-sound-speed anomalies if the sensing configuration only provides a lower ray coverage (i.e., < 100 ray per pixels for the studied sensing configurations). Finally, if the submesoscale variability is especially high (e.g., for sound speed anomalies > 1 m/s), such as for scenario C in this study, it is necessary to enhance the ray coverage to achieve good OAT estimates. In this case, adaptive methods can be beneficial, especially, if the additional acoustic transponders at each iteration can be strategically added to better probe the domain areas where the highest sound speed anomalies were presumed to be based on the previous iteration results. Hence, in practice, the ray coverage could be adjusted in an adaptive manner, potentially leveraging the maneuverability of small sensing platforms to best sample the expected submesoscale variability. Additionally, the constrained OAT estimator [Eq. (12)] was found to perform best for the adaptive OAT method, since it gradually yields a more distributed sets of constrained sound velocities locations across the inversion domain for increasing number of iterations. On the other hand, the unconstrained OAT estimator [Eq. (8)] yielded better estimates for the iterative OAT method. In this case, applying constraints only at a few discrete locations on the edges of the inversion domain (i.e., where the sparse acoustic transducers are located) created more estimation errors for the chosen formulation. In conclusion, the benefit of using a set of distributed sound velocity constraints (based on local measurements) to regularize the OAT inversion problem appears to be determined by the geometry of the actual sensing configuration. Figure 12 conceptually summarizes the applicability of the various OAT methods discussed in this study for different levels of environmental variability and acoustic sensor configurations.

Last, it is worth stressing that under the assumptions of this work, the estimation errors for the inverted SSMs represent only an idealized lower bound of what may be achievable by a practical OAT system with a similar acoustic sensing configuration. First and foremost, the submesoscale variability simulated by the ocean circulation model underestimates that of the real ocean due to horizontal and vertical resolution constraints. We plan to investigate this point in a future study. In addition, several sources of errors would likely increase the uncertainty of the OAT estimates in a practical situation, as extensively discussed in previous ray-based OAT literature,^7,14 including—among others—transducer positioning errors and unaccounted motion (e.g., array tilt due to tides or currents); imperfect time synchronization between the various source-receiver pairs and the impact of noise and available acoustic bandwidth—as well as the influence of ocean processes such as surface gravity waves, internal waves, and spice—on the accuracy of the measured ray travel-times.⁴⁶ Hence, the OAT performance is significantly driven by the available engineering resources which can ultimately enable the desired acoustic sensing configuration. However, the problem at stake in this study, namely, the estimation of the submesoscale ocean variability in the tactical upper ocean (< 1 km depth) over short ranges (< 10 km), may actually allow alleviating some of those engineering constraints. For example, for such short distances, the acoustic source sensors could operate at relatively high frequencies (> 10 kHz) with larger bandwidths, potentially providing a good signal-to-noise ratio and angle separation resolution which would ultimately improve the overall accuracy and diversity of the measured arrival times. Additionally, from an operational point of view, such compact high-frequency transducers could be more easily encapsulated in small (autonomous) profiling platforms, enabling simplified deployments and profiling of the upper ocean. These profiling platforms that would rely on rapid buoyancy change or other form or rapid mechanical actuation could provide depth-dependent sampling of the water column on small enough time scales to respect the frozen-ocean assumption between two OAT inversions. Furthermore, maintaining a surface expression (and hence GPS connection) for those autonomous profiling platforms may improve the power delivery to, and localization of, such profiling acoustic transducers, at least for the upper portion of the ocean domain of interest which typically exhibits the most dynamic submesoscale variability. In any event, sufficiently precise navigation systems should be implemented to ensure accurate positioning and tracking of the selected acoustic transponders, which is necessary for OAT applications.

This article compared the performance of different OAT methods over 5-day long estimation windows for three different scenarios exhibiting increasing submesoscale variability based on the simulated sound speed spatiotemporal variations computed using an ocean circulation model of the DeSoto Canyon in the Gulf of Mexico. Ray-based OAT were implemented for 2D upper ocean slices (depth < 1 km) and short ranges (< 5 km) for various acoustic sensing configurations that enabled various ray coverage density of the domain of interest. The results indicated that OAT methods have the potential to estimate submesoscale variability and highlighted the influence of the spatial and temporal scales of the sound speed anomalies of interest in the performance of the OAT estimation.

Additionally, the performance of iterative and adaptive OAT methods (encapsulating the classical OAT estimators) has been quantified to explore how to improve the estimation of sound speed anomalies in the upper ocean over short ranges when using a sparse or distributed sensor network, which could in practice be implemented with small sensing platforms. The results indicate that each of these OAT methods performs differently depending on the effective acoustic coverage and the sound speed variability to be estimated. Specifically, the classical OAT scheme was found to be sufficient for scenarios with high ray coverage and low sound speed variability, the iterated OAT scheme yielded better performances for lower acoustic coverage and low-to-moderate sound speed variability scenarios, and the adaptive OAT scheme was found to be most beneficial when the initial acoustic sensing configuration is sparse (thus yielding a low acoustic coverage initially), but the sound speed variability is high. These results could help guide the design of future OAT experiments that aim at validating the predictions from this numerical study and at quantifying the submesoscale variability in the upper ocean over short ranges.

The authors acknowledge the comments and suggestions made by the reviewers and the associate editor which significantly helped improve this manuscript. This work was supported by the Office of Naval Research (ONR) code 322 Ocean Acoustics Section, as well as the ONR Task Force Ocean (TFO) initiative under Grant No. N000142112558.

The authors have no conflict of interest to report.

Data supporting the findings of this study are available from the corresponding author upon reasonable request.

In practice, to implement the solution given by Eq. (A5), one possibility is to simply set $c const = c true$⁠, where $c true$ is the true SSM to be estimated and also corresponds to the sound speed values, which can be independently measured in situ at the selected locations N_c. However, it should be noted that in this study, the reference SSM $c 0$ [set to be the last SSM of the training set, just before the estimation period starts (see Fig. 2) used to compute the linearized forward model E [see Eqs. (1)–(4)] actually differs from the averaged SSM $c Φ ¯$ of the training set used to compute EOF basis. Hence, an alternate choice for the constrained values is $c const = c true − ( c Φ ¯ − c 0 )$ as implemented in this study in Eq. (12). The second term $( c Φ ¯ − c 0 )$ allows us in a simple fashion to account for the evolution of the background SSM that has occurred between $c Φ ¯$ and $c 0$ due to the local ocean dynamics. In this case, Eq. (A5) becomes identical to Eq. (12) which was used as the multi-objective OAT estimator in Sec. IV of this study such that the last penalty term effectively reduces to $λ 2 U Φ T U ( c true − c 0 ) ]$ (which is referred to as the $c 0$ case) instead of $λ 2 U Φ T U ( c true − c Φ ¯ ) ]$ (which is referred to as the $c Φ ¯$ case).

Figure 13 compares the performance of these two types of constrained OAT for scenarios A and C based on the same adaptive setup algorithm used throughout this study and presented in Fig. 5(b). Overall, these results indicate the formulation for the multi-objective OAT estimators given by Eq. (12) [i.e., by setting $c const = c true − ( c Φ ¯ − c 0 )$ in the general formulation of Eq. (A5)]—corresponding to the dashed-dotted lines—performs slightly better for the selected test parameters and sensing configurations.