{'fraction_non_alphanumeric': 0.048499876022811805, 'fraction_numerical': 0.02355566575750062, 'mean_word_length': 4.384779706275033, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimating ocean's density stratification from surface data

Here we propose a semi-analytical technique that accurately reconstructs ocean's density stratification profile, and hence, the pycnocline depth, simply from the free surface elevation data. The ocean surface contains the signature of internal gravity waves (IGWs), which are generated when stably stratified ocean water is forced to move back and forth over submarine topography by barotropic tides. Since IGWs, in turn, contain the information of ocean's density stratification, the latter can in principle be reconstructed from the free surface signature. We consider a simple theoretical model that approximates a continuously stratified ocean as discrete layers of constant buoyancy frequency; this facilitates the derivation of a closed-form dispersion relation. First, we numerically simulate IGW generation for toy ocean scenarios and subsequently perform space-time Fourier transform (STFT) of the free surface. Free surface STFT yields IGW spectra which has wavenumbers corresponding to the tidal frequency. Density profile reconstruction has been performed by substituting these IGW wavenumbers into the dispersion relation. Next we consider a more realistic situation in which rotation is included; moreover the bottom topography, density and shear profiles are representative of the Strait of Gibraltar. While density reconstruction is 94.4% accurate in the absence of shear, the accuracy decreases to 90.2% when moderate shear is present. Since shear causes loss of coherence in IGW beams, density reconstruction can be adversely affected in those regions of the global ocean where shear is very strong. Finally, using wavenumbers obtained from multi-satellite altimetry data near the Hawaiian Ridge, we have reconstructed actual density profile with 94.4% accuracy. Since shear is either low or moderate in most parts of the global ocean, we expect our proposed technique to be broadly applicable.


Introduction
Oceans are by and large stably stratified, that is, the density of ocean water monotonically increases with depth. Ocean's density also varies with latitude and longitude, as well as with seasons. Depending upon the strength of stratification, the vertical structure of ocean's density is divided into three major layers: (i) top -weakly stratified surface mixed layer, (ii) middle -strongly stratified pycnocline, and (iii) bottom -weakly stratified abyss (Sutherland 2010).
An accurate knowledge of ocean's density field is crucial for ocean and climate modeling (Cummins 1991). One of the most important consequences of stable density stratification is internal gravity waves (IGW). In oceans, IGWs are often produced when the stably stratified ocean water is driven back and forth over submarine topography by tidal currents; such IGWs are also known as internal tides. In general, IGWs are very efficient in transporting momentum and energy over large distances, and help in mixing of nutrients, oxygen and heat in the oceans (Sarkar & Scotti 2017). Turbulence and mixing via IGWs play an important role in regulating global oceanic circulation, and are one of the major factors in the climate-forecast models (Ferrari & Wunsch 2009). Oceanic density stratification also has a direct impact on the aquatic ecosystem. In oceans and lakes, microbiological activities and accumulation of organisms are strongly affected by the pycnocline (Doostmohammadi et al. 2012). Density stratification influences the formation of phytoplankton blooms, which in turns help to maintain a balanced ecosystem (Sherman et al. 1998).
Ocean's density is a function of both temperature and salinity, both of which are measured using CTD (Conductivity, Temperature and Depth) sensors. These sensors, while descending (or ascending) through the ocean water, collects the necessary information. The vertical profiles of temperature and salinity thus obtained are then substituted into the equation of state to yield ocean's density profile at a given latitude-longitude. To the best of our knowledge, there are no indirect or non-invasive techniques that can estimate oceanic density profile.
In this article, we propose a strategy that can provide a reasonably accurate estimate of ocean's density stratification profile (and hence, the pycnocline depth) in a fully noninvasive manner by only analyzing the ocean free surface. Satellite images have already revealed that in regions with strong tidal forcing, the free surface gets perturbed by internal tides (Ray & Mitchum 1996). These perturbations travel as surface waves and their frequency matches the tidal frequency. Recently, efforts have been made to study these waves using in-situ measurements and satellite altimetry, in conjunction with global circulation models (Ray & Zaron 2011;Savage et al. 2017a,b). The pivotal point of this article is the realization that these surface signatures carry the information of ocean's density stratification, and can in principle be inverted to reconstruct the latter.
The paper is organized as follows. The theoretical and numerical aspects are respectively discussed in §2 and §3. In §4, we first consider toy models with simple density profiles, and finally, a representative density profile of the Mediterranean sea. In each case, IGWs emanating from the bottom topography impinges on the free surface. Wavenumbers corresponding to the surface signature are substituted in a closed form dispersion relation to reconstruct the underlying density profile. The article is concluded in §5.

Governing equations and exact solutions
We consider an incompressible, 2D (x−z plane), density stratified flow of a Boussinesq fluid. The mean (denoted by overbars) density profile varies in the vertical (z) direction. The governing Navier-Stokes equation in this case is given by (Gerkema 2001;Gerkema & Zimmerman 2008;Verma 2018): Except F, variables without overbars denote perturbation quantities. The perturbation velocity field is denoted by u ≡ (u, w), p and ρ respectively denote perturbation pressure and density. The quantity g denotes gravitational acceleration, ρ 0 represents reference density, ν is the kinematic viscosity and κ is the mass diffusivity. Furthermore, N (z) ≡ −(g/ρ 0 )dρ/dz is the Brunt-Väisälä (or buoyancy) frequency, which is a measure of the background stratification. The forcing function, F, is assumed to be zero here. In the linear regime, (2.1a)-(2.1c) can be simplified into one equation by neglecting the effect of viscosity and diffusivity, and can be expressed as We seek a plane-wave solution and express w as follows: w = W (z)e i(kx−ωt) , where k is the wavenumber in the x-direction and ω is the frequency. By substituting this ansatz in (2.2), we get We assume the lower boundary (at z = −H) to be impenetrable, i.e. w = 0 (implying W = 0), while the upper boundary (z = 0) to be a free surface. However, as shown in Appendix A, the leading order approximation yields w = 0 (implying W = 0) even at the free surface -which is popularly known as the 'rigid-lid approximation'. Equation (2.3), together with the homogeneous boundary conditions, constitute a regular Sturm-Liouville boundary value problem. Its solution is formed by the superposition of a countably infinite set of eigenvalues k n and corresponding eigenfunctions W n . The solution, which physically represents internal gravity waves, can be obtained in analytical form only for some special choices of N (z), e.g. constant or piecewise constant (Gerkema & Zimmerman 2008), otherwise (2.3) has to be solved numerically. Below we provide the exact solutions for (i) a single layer of constant N , (ii) two layers, each having a constant N , and (iii) three layers, each having a constant N .

One layer
We consider a mean density profileρ(z) that varies linearly with z, giving a constant N . In this situation, (2.3), along with the homogeneous Dirichlet boundary condition, can be solved exactly, yielding where C n ∈ R are arbitrary constants. The dispersion relation is given by (2.5) In the above equation, we have fixed the value of ω as ω 0 , which we take as the tidal frequency. This is because our interest here is to obtain internal tides, that is, IGWs oscillating at tidal frequencies. Equations (2.4) and/or (2.5) appear in classic texts, e.g. Turner (1979); Gerkema & Zimmerman (2008); Sutherland (2010).

Two layers
In this case we consider a two-layered density stratified flow. Density in each layer varies linearly with z (i.e., N is constant in each layer) as follows: We note here that the two-layered density stratification is such thatρ(z) is still continuous, implying there are no interfacial gravity waves at the pycnocline z = −h. In such a system, W n can be written as (Gerkema & Zimmerman 2008): By demanding the continuity of W n and dW n /dz at z = −h, we arrive at the dispersion relation where m n,i = k n N 2 i − ω 2 0 /ω 0 ; i = 1, 2.

Three layers
Next we consider a three-layered density stratified flow, with N = constant in each layer: (2.8) Again we note that theρ(z) is continuous. As already mentioned, oceans can be broadly divided into three regions of different density stratifications, hence the three-layered model can crudely represent ocean's mean density profile. Thus N 1 , N 2 and N 3 respectively denote the stratifications of the top, middle (pycnocline) and bottom layers. In this three-layered system, W n can be expressed as where m n,i = k n N 2 i − ω 2 0 /ω 0 ; i = 1, 2, 3. To obtain the four unknown coefficients, we demand the continuity of W n and dW n /dz at the two interfaces z = −h 1 and z = −h 2 , which finally yields the dispersion relation

Numerical implementation
In order to simulate internal tides, we numerically solve (2.1a)-(2.1c). Following Gerkema (2001), we consider the barotropic tidal forcing term where Q is the flow rate, h(x) is the local water depth and ω 0 is the tidal frequency. Bottom topography has been incorporated; furthermore, the Cartesian coordinate system, x − z, has been transformed into a terrain-following coordinate system, x − ζ with ζ ≡ −z/h(x). Therefore, in the terrain-following coordinate system, the undisturbed free surface is denoted by ζ = 0, and the bottom surface lies at ζ = −1. Following Dimas & Triantafyllou (1994), we use a spectral spatial discretization with Chebyshev polynomial in the vertical direction and Fourier modes in the streamwise direction, the latter has been assumed to be periodic. Equations (2.1a)-(2.1c) have been solved using an opensource pseudo-spectral code -Dedalus (Burns et al. 2017). Fourth-order Runge-Kutta method has been used for time-marching.
The correspondence between the prognostic variables in the two coordinate systems are:ũ(x, ζ, t) = u(x, z, t);p(x, ζ, t) = p(x, z, t);ρ(x, ζ, t) = ρ(x, z, t). At the free surface ζ = η(x, t), the kinematic and dynamic boundary conditions are respectively given bỹ where η is the free surface elevation. At the free surface,ũ satisfies the stress free boundary condition. At the bottom surface ζ = −1,ũ andw respectively satisfies the no-slip and no-penetration boundary conditions. Insulating boundary conditions both at top and bottom have been used for density. Furthermore, we have taken both viscosity and diffusivity into account, ν and κ are respectively set to 10 −6 m 2 s −1 and 10 −7 m 2 s −1 .
For numerical simulations, we have first considered a toy model with a Gaussian bottom topography. The density profile of the model has been varied from single to three-layered (given in §4.1- §4.3). The model has a depth of H = 1 m, a horizontal extent of 20 m, and has been forced with a barotropic tidal flow of amplitude 10 −3 m s −1 and frequency ω 0 = 0.05 s −1 . Since the topography radiates IGWs, sponge layers have been used to absorb the incoming IGWs both at the east and the west boundaries of the domain. We have used 256 Chebyshev points in the z-direction and 1024 Fourier-modes along the x-direction. We have simulated 8 tidal periods with a time-step of 0.1 s.
Next we have considered a realistic scenario in §4.4 and studied IGW generation in the tidally active part of the Mediterranean sea. The stably stratified, time-averaged and smoothed density profile has been taken at 36.6 • N latitude and 0.2 • W longitude. To simulate this we have used a domain of (L x × L z ) = (50 × 1) km with sponge layers of 10 km on both eastern and western boundaries. The flow is forced using semi-diurnal tides of frequency 1.4 × 10 −4 s −1 and amplitude 10 −3 ms −1 over a Gaussian mountain. The spatial and temporal discretization, as well as the total time are same as that in §4.1- §4.3.

One-layer
We first simulate a single-layered flow with N = 0.1 s −1 . Our objective is to get back this value of N by only analyzing the free surface data. Due to the tidal forcing, internal wave beams radiate from the Gaussian topography and impinges on the free surface. The space-time Fourier Transform (STFT) of the surface elevation field η yields the wavenumbers k n corresponding to the tidal frequency ω 0 = 0.05 s −1 . Figure 1(a) shows STFT of the free surface elevation; the first vertical-mode (n = 1) corresponds to k 1 = 1.813 m −1 . By substituting ω 0 and k 1 in (2.5), we straightforwardly estimate N . Since the density at the surface is known, the mean density profileρ(z) can be directly reconstructed; see figure 1(b). Moreover, this result also serves as a validation of the numerical code.

Two-layers
Here we consider a squared buoyancy frequency N 2 = 2 × 10 −2 − 10 −2 1 1 + (z − 0.7) 256 , (4.1) which closely resembles (2.6) with N 1 = 0.1 s −1 , N 2 = 0.14s −1 and h = 0.3 m. Similar to the one-layer case in §4.1, tidal forcing leads to IGW radiation, whose imprint is detectable at the free surface. Figure 1(c) shows STFT of the free surface elevation; corresponding to the tidal frequency, the first three vertical modes respectively peak at k 1 = 1.813 m −1 , k 2 = 2.549 m −1 and k 3 = 3.836 m −1 . Our objective is to reconstruct (4.1) using (2.7), which means that the three unknowns, N 1 , N 2 and h have to be evaluated. We substitute the obtained values of k 1 , k 2 and k 3 along with ω 0 in (2.7), leading to a system of three equations and three unknowns, which is then solved numerically. The mean density profile along with the reconstructed version are shown in figure 1(d).

Three-layers
We follow the same strategy as that outlined in §4.1 and §4.2. In this case, we use the following N 2 profile: (4.2) The above profile closely resembles the three-layered configuration (2.8), with N 1 = 0.1 s −1 , N 2 = 0.25s −1 , N 3 = 0.1s −1 , h 1 = 0.12 m and h 2 = 0.28 m. The goal is to find N 1 , N 2 , N 3 , h 1 and h 2 , and therefore we construct five-equations from (2.9) for different value of n. Figure 1(e) reveals that k 1 = 1.596 m −1 , k 2 = 2.448 m −1 , k 3 = 4.789 m −1 , k 4 = 6.173 m −1 and k 5 = 8.408 m −1 . The system of five equations and five unknowns resulting from (2.9) are then solved numerically. The actual and the reconstructed profiles are shown in figure 1(f).

Mediterranean sea profile
As already mentioned in §3, we simulate IGWs for a case in which the mean density profile is representative of the Mediterranean sea. Figure 2 (also see supplementary Movie 1) shows snapshots of the free surface displacement η along with the contours Figure 2: IGWs radiating from a Gaussian submarine mountain impinges on the free surface. Density profile, representative of the Mediterranean sea, has been considered. Snapshots of the free surface displacement (η, in m) and the corresponding horizontal baroclinic velocity field (u, in ms −1 ) are shown. Time corresponding to each snapshot appears at the top of each sub-figure. The semi-diurnal tidal period, T = 2π/ω ≈ 12.46 hours. Bending of the internal beams occur due to refraction from the pycnocline, furthermore, reflection from the pycnocline (which is of moderate strength) leads to the observed beam scattering. A point worth mentioning here is that the wavenumbers present at the free surface displacement (and hence, in the internal beam) do not depend on the underlying topography at the generation site. Therefore, the result should be valid for any other bottom topography that is not flat. Figure 3(a) represents STFT of the free surface elevation, the first five vertical-modes are k 1 = 1.49 × 10 −4 m −1 , k 2 = 3.19 × 10 −4 m −1 , k 3 = 5.32 × 10 −4 m −1 , k 4 = 6.81 × 10 −4 m −1 and k 5 = 8.94 × 10 −4 m −1 . We approximate the Mediterranean sea profile with a three-layered model, and hence follow the same procedure mentioned in §4.3. The reconstructed density profile shown in figure  3(b) has 99.9% accuracy, implying that the three-layered model (or in other words, the first five modes) is sufficient in this regard.

Conclusion
In this paper we numerically show that ocean's free surface carries information regarding its mean density profile, and therefore devise a semi-analytical strategy towards its reconstruction. To the best of our knowledge, this is the first work that puts forward an invasive technique towards obtaining or estimating ocean's density profile. Barotropic tides cause stratified ocean water to move back and forth over submarine topography, causing IGW radiation. IGWs carrying information regarding ocean's mean density profile impinge on the free surface, the signature of which can be observed using STFT. Wavenumbers (countably infinite in number) constituting this IGW correspond to the tidal frequency in the STFT spectrum. In general, higher the mode number, lower is its amplitude; hence the first few modes of the surface signature (which will be easier to detect in practice) can be used to estimate the density stratification profile.
Ocean's mean density profileρ(z) has some specific qualities -it continuously and monotonically decreases with z. Furthermore, the variation ofρ(z) is such that it can be broadly divided into three distinct regions. In each of these regions,ρ(z) can be approximated to be linearly varying with z (implying N is constant in each layer). Therefore we construct a 'simplified ocean' with m layers, each with a constant N (the maximum value of m considered here is 3). The remarkable advantage of this simplification is that the a closed form dispersion relation can be obtained. For an m-layered flow (m 2) we have to find buoyancy frequencies in each layer: N 1 , N 2 , · · · , N m , and layer depths: h 1 , h 2 , · · · , h m−1 . It is possible to evaluate these 2m − 1 unknowns by constructing 2m − 1 equations out of the dispersion relation, provided we know the wavenumbers k 1 , k 2 , · · · , k 2m−1 corresponding to the tidal frequency ω 0 . Indeed, this information is known from the STFT spectrum of the free surface. After reconstructing simpler density profiles, we consider a more complicated density profile that is representative of the Mediterranean sea. Using a 3-layered model, we reconstruct the above-mentioned profile with 99.9% accuracy. The reconstruction accuracy is expected to be lower in real ocean scenarios since the signature of the higher modes at the free surface may not be detectable. In fact, the first step towards actual density profile reconstruction would be to capture the first five modes, which would yield the three-layered structure. The upcoming satellite mission Surface Water and Ocean Topography (SWOT) may have the capabilities to detect these high wavenumber sea surface height signatures (Savage et al. 2017a).
velocity scales: c, which is the phase speed of internal gravity waves, and c sg , which is the long surface gravity wave speed. Therefore, the non-dimensional quantity shows ratio between phase speed of the internal gravity waves and large-wavelength surface gravity waves. Equation (2.5) reveals the characteristic c is ≈ N H/nπ (when N ω 0 ). For a typical ocean of depth 2 km, taking N = 10 −3 s −1 and n = 1 gives c ≈ 0.63 ms −1 and c sg ≈ 140 ms −1 . Therefore ≈ 2 × 10 −5 O(1). Hence at the leading order we have w = 0, (A 8) which means that the free surface at the leading order acts like a rigid-lid boundary.