Probabilistic Ocean Acoustic Modeling and Coupled Bayesian Data

Probabilistic Ocean Acoustic Modeling and Coupled Bayesian Data Assimilation for Bathymetry Inversion M. S. Bhabra, W. H. Ali and P. F. J. Lermusiaux January 2020 Multidisciplinary Simulation, Estimation, and Assimilation Systems (MSEAS) Department of Mechanical Engineering, MIT http: //mseas. mit. edu 1

Main MIT‐MSEAS contributions (August 2019 ‐ January 2020) • Completed and submitted a joint paper for IEEE OCEANS 2019. • Further developed the acoustic wavefront codes to model scattering off arbitrary curved boundaries. • Derived wavefront propagation equations that efficiently account for uncertainties in the sound speed field and seabed bathymetry. • Developed and implemented a ray‐tracing code including its preliminary DO‐ stochastics version. • Made good progress on the bathymetry inversion framework and implemented techniques information‐based techniques to assist in locating the receivers for better performance. • Set up ocean and acoustic modeling for the Mass Bay area in preparation for sea tests. 2

Background: Underwater Acoustic Modeling • Homogeneous Acoustic Wave Equation (Godin, Springer‐ 1987): Time‐domain methods Frequency‐domain methods ‐ Finite Difference, Finite Element methods ‐ Level‐set methods ‐ ‐ Ray methods Wavenumber integration Normal Modes Parabolic Equation 3

Background: Underwater Acoustic Modeling • Homogeneous Acoustic Wave Equation (Godin, Springer‐ 1987): Time‐domain methods Frequency‐domain methods ‐ Finite Difference, Finite Element methods ‐ Level‐set methods ‐ ‐ Ray methods Wavenumber integration Normal Modes Parabolic Equation 4

The Level Set Method for Wavefront Propagation • Consider high frequency solutions to the Acoustic Wave Equation using a solution of the following form: • is the phase function. At any given point of time, the level surfaces of correspond to points at which the phase is the same (surfaces of constant phase, also referred to as the wavefront). • Upon substituting into the wave equation (and retaining the highest order terms in ) obtain (Martinelli, 2012): • Has the form of a Hamilton‐Jacobi Equation 5

The Level Set Method for Wavefront Propagation (2) • Two possible approaches can be used to solve the Hamilton‐Jacobi equation dictating the wavefront evolution Lagrangian Approach • Considers a discrete set of particles on the initial wavefront: • Trajectory of each point governed by characteristics of the Hamiltonian system (Engquist et al. , 2002): Eulerian Approach • Considers the front to be the level surface of an implicit function in a higher dimension: • Evolve the implicit function by solving a PDE governing the function’s evolution (the Hamilton‐Jacobi equations). • Known as the Level Set Method • Is a Ray Tracing approach • Difficulty: Points originally close together may diverge at later times (poor wavefront resolution) • To evolve the wavefront, codimension‐two objects are used (to allow for multivalued solutions): • Represent the front in the level set framework as the intersection of the isocontours of two functions. 6

The Level Set Method for Wavefront Propagation (3) • Evolve the acoustic wavefront as a curve/strip in a higher dimensional reduced phase space with coordinates (�� , �� ). • �� gives the normal direction of a given point on the wavefront (measured from the +�� axis). • The wavefront is found by projecting this curve onto the (�� , �� ) plane. • The velocity field for the strip’s propagation can be derived to be (Engquist et al. , 2002): • The strip is represented as the intersection of two level set functions whose evolution is governed by (Martinelli, 2012; Osher et al. , 2002) : • The wavefront (surfaces of constant phase) is given as: 7

The Level Set Method for Wavefront Propagation (4) • To handle non‐rectangular physical domains, a mapped form of the physical equations is solved along with proper treatment of the boundary for efficient specular reflection. • To efficiently and accurately parametrize the physical domain, a spline‐based approach has been used in which the domain is transformed using a set of control points and basis functions to ultimately yield a mapping of the form: • Due to the cartesian nature of the phase dimension in the augmented space, no special mapping is needed here (an identity mapping may be used). • The final mapping from the augmented parametric domain to the physical domain is then given as: 8

The Level Set Method for Wavefront Propagation (5) • The transformation may then be used to obtain the mapped equations governing the level set evolution in the parametric domain: where is the determinant of the Jacobian for the mapping . • The specular reflection boundary conditions is moreover implemented by using the grid transformation as: • The mapped equations allow for an efficient means to handle arbitrary domains. • Discretization Details: • Temporal Terms: Can discretize explicitly in time. In this work, have used TVD Runge‐Kutta Schemes. • Spatial Terms: Discretize in space using high order Essentially Non‐Oscillatory (ENO) or Weighted Essentially Non‐Oscillatory (WENO) operators 9

Case 1: Variable Wave Speed (Linear Boundary) Spatial Discretization Temporal Discretization Time Stepping Scheme Second Order TVD Runge‐Kutta Spatial Derivatives Scheme Second Order ENO Reflecting Boundaries 10

Case 2: Variable Wave Speed (Curved Boundary) Spatial Discretization Temporal Discretization Time Stepping Scheme Third Order TVD Runge‐Kutta Spatial Derivatives Scheme Third Order ENO Reflecting Boundaries 11

Case 3: 12

Parabolic Equation (PE) Modeling HE PE Helmholtz Equation Modeling Parabolic Equation Modeling 13

Parabolic Equation (PE) Modeling • “Standard” Narrow-Angle PE (NAPE) (Tappert, Springer‐ 1977): • T‐s approximation of 1 st order with propagation angle • Obtain reaction‐diffusion like equation, • Wide-Angle PE (WAPE) (Sturm and Fawcett, JASA‐ 2003): • P‐s approximation Padé(m): • Substitute in PE: • Padé(1) corresponds to the Claerbout WAPE with propagation angle Implemented range-dependent NAPE & WAPE (2 nd order FV & range-marching) • New efficient Padé(m) WAPE that solves m intermediate prob. (instead of inverting operator) • Verified using classic test cases & comparisons (KRAKEN, BELLHOP, etc. ) in realistic oceans 14

Validation for 2 D Range‐Independent and Range‐Dependent Benchmark Problems [Davis et al. , NORDA‐ 1982; Jensen and Ferla, JASA‐ 1990] WAPE 15

Evaluation for Realistic Ocean Test Cases • Ocean sound speed section: • • Across the shelf break in the Middle Atlantic – New York Bight Region Obtained from a forecast by the data assimilative MSEAS primitive equation model during a DARPA POSYDON‐POINT real‐time sea experiment in August, 2018 Sound source located within the sound channel Using KRAKEN was computationally infeasible due to complex range dependency BELLHOP (Porter, NRL) 16

Stochastic Dynamically‐Orthogonal NAPE Modeling • Stochastic NAPE: • Substitute DO representation [Sapsis and Lermusiaux, 2009, 2011; Ueckermann et al. , 2013]: similarly Mean PDE Modes PDEs ⇒ Coefficients SDEs [Ali and Lermusiaux, 2020 a] 17

Bayesian Ocean Physics‐Acoustics‐Bathymetry DA – Beyond Tomography State Vector Augmentation Approach • Augment ocean physics variables with acoustics variables • Augmented State: § Augmented dynamical model: § Diagnostic variable: v Transmission Loss: Bayesian DA Algorithm: Nonlinear smoothing – all ranges at once • Full Bayesian inverse: using GMM‐DO Filter [Sondergaard and Lermusiaux, MWR‐ 2013 a, b; Lolla and Lermusiaux, MWR‐ 2017 a, b] § § Capture non‐Gaussian distributions Perform assimilation in the reduced subspace from DO Assimilate measurements at particular ranges Correct predictions at all ranges 18 [Ali and Lermusiaux, 2020 b]

2 D Bump Propagation (Ali et al. , 2019) 26

Evolution of TL Mean, Std. Dev. and Modes Number of modes = 30, number of MC samples = 1000 27

Bayesian GMM‐DO DA of Sparse TL Data Objective: Jointly learn range‐dependent sound‐speed, bathymetry, and TL fields by assimilating sparse TL data TL Measurements (10 Total) Difference

Bayesian GMM‐DO DA of Sparse TL Data Objective: Jointly learn range‐dependent sound‐speed, bathymetry, and TL fields by assimilating sparse TL data TL Measurements (10 Total) Difference

Bayesian GMM‐DO DA of Sparse TL Data Objective: Jointly learn range‐dependent sound‐speed, bathymetry, and TL fields by assimilating sparse TL data

OCEANS 2019 paper 26

Ongoing Work 25

References Ali, W. H. and P. F. J. Lermusiaux, 2019 a. Dynamically Orthogonal Equations for Stochastic Underwater Sound Propagation: Theory, Schemes and Applications. In prep. Ali, W. H. and P. F. J. Lermusiaux, 2019 b. Generalized Ocean Acoustics Bayesian Inversion with Gaussian Mixture Models using the Dynamically Orthogonal Field Equations. In prep. Brooke, G. H. , D. J. Thomson, and G. R. Ebbeson, 2001. PECan: A Canadian parabolic equation model for underwater sound propagation. Journal of Computational Acoustics, 9(01), pp. 69‐ 100. Chin‐Bing, S. A. , D. B. King, and J. E. Murphy, 1993. Numerical simulations of lower‐frequency acoustic propagation and backscatter from solitary internal waves in a shallow water environment. In Ocean Reverberation (pp. 113‐ 118). Springer, Dordrecht. Claerbout, J. F. , 1985. Fundamentals of geophysical data processing with applications to petroleum prospecting. Mc. Graw‐Hill, New York. Davis, J. A. , D. White, and R. C. Cavanagh, 1982. Norda parabolic equation workshop, 31 March‐ 3 April 1981 (No. NORDA‐TN‐ 143). Naval Ocean Research and Development Activity Stennis Space Center MS. Engquist, B. , O. Runborg, and A. K. Tornberg, 2002. High‐frequency wave propagation by the segment projection method. Journal of Computational Physics, 178(2), pp. 373‐ 390. Jensen, F. B. and C. M. Ferla, 1990. Numerical solutions of range‐dependent benchmark problems in ocean acoustics. The Journal of the Acoustical Society of America, 87(4), pp. 1499‐ 1510. Lee, D. and S. T. Mc. Daniel, 1988. Ocean acoustic propagation by finite difference methods (Vol. 15). Elsevier. 26

References Lin, Y. T. , Duda, T. F. and A. E. Newhall, 2013. Three‐dimensional sound propagation models using the parabolic‐equation approximation and the split‐step Fourier method. Journal of Computational Acoustics, 21(01), p. 1250018. Lolla, T. and P. F. J. Lermusiaux, 2017 a. A Gaussian mixture model smoother for continuous nonlinear stochastic dynamical systems: Theory and scheme. Monthly Weather Review, 145(7), pp. 2743‐ 2761. Lolla, T. and P. F. J. Lermusiaux, 2017 b. A Gaussian mixture model smoother for continuous nonlinear stochastic dynamical systems: Applications. Monthly Weather Review, 145(7), pp. 2763‐ 2790. Martinelli, S. L. , 2012. An application of the level set method to underwater acoustic propagation. Communications in Computational Physics, 12(5), pp. 1359‐ 1391. Osher, S. , L. T. Cheng, M. Kang, H. Shim, and Y. H. Tsai, 2002. Geometric optics in a phase‐space‐based level set and Eulerian framework. Journal of Computational Physics, 179(2), pp. 622‐ 648. Sapsis, T. P. and P. F. J. Lermusiaux, 2009. Dynamically orthogonal field equations for continuous stochastic dynamical systems. Physica D: Nonlinear Phenomena, 238(23‐ 24), pp. 2347‐ 2360. Sapsis, T. P. and P. F. J. Lermusiaux, 2012. Dynamical criteria for the evolution of the stochastic dimensionality in flows with uncertainty. Physica D: Nonlinear Phenomena, 241(1), pp. 60‐ 76. Subramani, D. N. , Q. J. Wei, and P. F. J. Lermusiaux, 2018. Stochastic Time‐Optimal Path‐Planning in Uncertain, Strong, and Dynamic Flows. Computer Methods in Applied Mechanics and Engineering, 333, 218– 237. doi: 10. 1016/j. cma. 2018. 01. 004 27

References Tappert, F. D. , 1977. The parabolic approximation method. In Wave propagation and underwater acoustics (pp. 224‐ 287). Springer, Berlin, Heidelberg. Ueckermann, M. P. , P. F. J. Lermusiaux, and T. P. Sapsis, 2013. Numerical schemes for dynamically orthogonal equations of stochastic fluid and ocean flows. Journal of Computational Physics, 233, pp. 272‐ 294. 28