Joint Estimation of Balanced Motions and Internal Tides From Future Wide-Swath Altimetry

a


Introduction
For more than 25 years, altimetry has allowed the study of near-global sea surface height (SSH) at scales longer than 150 km and drastically transformed our understanding of mesoscale processes in the oceans. In particular, altimetry has given access to the upper-ocean circulation dynamics through geostrophy and revealed that most kinetic energy in the world ocean is contained in mesoscale eddies with wavelengths between 100 and 300 km (Wunsch & Ferrari, 2004). But the scales resolved by nadir (along-track) altimetry are limited by the 100-300 km spacing between satellite ground tracks. Even by merging several nadir measurements , the space time resolution of the resulting 2D SSH maps does not allow to conveniently characterize and study small mesoscale (<150 km) motions (Amores et al., 2018;Ballarotta et al., 2019).
With wide-swath radar interferometry, the next generation of altimeters will provide two dimensional SSH data at much finer scales. For instance, the future Surface Water and Ocean Topography (SWOT) mission (Morrow et al., 2019) will provide SSH observations over a 120 km-wide swath with a near-global coverage, resolving spatial scales down to 15-30 km depending on the sea state (Fu & Ubelmann, 2014). Recent studies, based Abstract Wide-swath altimetry, for example, the Surface Water and Ocean Topography mission is expected to provide Sea Surface Height (SSH) measurements resolving scales of a few tens of kilometers. Over a large fraction of the globe, the SSH signal at these scales is essentially a superposition of a component due to balanced motions (BMs) and another component due to internal tides (ITs). Several oceanographic applications require the separation of these components and their mapping on regular grids. For that purpose, the paper introduces an alternating minimization algorithm that iteratively implements two data assimilation techniques, each specific to the mapping of one component: a quasi-geostrophic model with Back-and-Forth Nudging for BMs, and a linear shallow-water model with 4-Dimensional Variational assimilation for ITs. The algorithm is tested with Observation System Simulation Experiments where the truth is provided by a primitive-equation ocean model in an idealized configuration simulating a turbulent jet and mode-one ITs. The algorithm reconstructs almost 80% of the variance of BMs and ITs, the remaining 20% being mostly due to dynamics that cannot be described by the simple models used. Importantly, in addition to the reconstruction of stationary ITs, the amplitude and phase of nonstationary ITs are reconstructed. Sensitivity experiments show that the quality of reconstruction significantly depends upon the timing of observations. Although idealized, this study represents a step forward towards the disentanglement of BMs and ITs signals from real wide-swath altimetry data.
• We present a dynamical method that estimates and separates balanced motions (BMs) and internal tides (ITs) from sea surface height data • The method uses iteratively two data assimilation techniques, each specific to the estimation of one component (BMs or ITs) • Although idealized, this study shows encouraging results for the disentanglement of BMs and ITs signature on future wide-swath altimetry data 2 of 17 on numerical models, have highlighted the impacts of short-mesoscale and submesoscale processes on ocean dynamics. Submesoscale motions have been found to trigger a large part of vertical motions (Klein & Lapeyre, 2009;McWilliams, 2016) that drive the exchanges of heat, carbon and nutrients between the ocean surface and subsurface. Small mesoscale processes (<100 km) also provide kinetic energy to larger scales through an inverse cascade (Ajayi et al., 2019;Capet et al., 2016), hence impacting the mesoscale dynamics. Wide-swath altimetry represents a unique opportunity to validate our present-day understanding of submesoscale motions and their role in oceanic processes in a broad sense.
At the scales measured by wide-swath altimetry, some SSH variations due to internal tides (ITs) may become comparable to those due to the mesoscale balanced motions (BMs; Qiu et al., 2018). ITs, producing high frequency sea level fluctuations at scales around 200 km and below, are internal gravity waves generated when the barotropic tidal flow encounters variations of topography (Garrett & Kunze, 2007). Although the generation processes of ITs are well understood, the dissipation processes (that greatly influence the ocean's energy budget; e.g., Whalen et al., 2020) remain insufficiently known and quantified. Recent studies emphasize that ITs and BMs interact in a complex way (Kelly & Lermusiaux, 2016), with BMs scattering and refracting ITs . These complex interactions are thought to strongly modify local mixing but need to be validated with observations. The capacity of wide-swath altimetry to observe even a part of ITs SSH signal will help to better understand these interactions and the processes related to ocean mixing and dissipation.
The accurate characterization of both ITs and BMs from wide-swath altimetry data, in the small mesoscale spectrum where their signatures overlap, will require scientific and technological developments for the processing of observations. The first challenge concerns the design of gridded products. For BMs, the main objective is to increase the space-time resolution of the present-day SSH maps. For instance, the mismatch between the (high) spatial and (low) temporal sampling of SWOT fosters the development of innovative inversion techniques, which typically include dynamical constraints (Ubelmann et al., 2015). For ITs, the main objective is to map the incoherent part of the signal, that is, the part due to waves whose characteristics have been altered by interactions with BMs and stratification variations. This is further developed in the next paragraph. The second challenge consists in the separation of ITs and BMs components from wide-swath altimetry data to make the differentiated gridded products possible. Both components are driven by different dynamics, and separation is needed to properly estimate quantities associated with each (e.g., surface currents for tracer advection, or energy dissipation due to waves). Separation is also made difficult by geographical and seasonal variations in the spatial scales and the relative strength of ITs and BMs (Qiu et al., 2018).
Disentangling the contribution of BMs and ITs on SSH is still an unsolved challenge. Due to its partial sampling in space and time, altimetry does not capture the fast SSH pulsation due to ITs. Stationary ITs (which are phaselocked to the astronomical forcing) can be mapped by harmonic analysis of long time series (>10 years) of conventional nadir altimetry data (Ray & Zaron, 2016). ITs become non stationary when interacting with BMs Ponte & Klein, 2015) and/or being modulated by the stratification (Ray & Zaron, 2011). The phase and amplitude modulations of the ITs field can evolve on time scales between 5 and 20 days, which makes their predictability nearly impossible with conventional satellite altimetry (Haren et al., 2004;Nash et al., 2012). Recent studies (Nelson et al., 2019;Zaron, 2017) suggest that nonstationarity represents half of the total ITs variance on average. Wide-swath altimetry will considerably increase the SSH measurement density and open the way to predicting nonstationary ITs.
In the context of SWOT development, the BMs-ITs separation problem has been addressed with various approaches. Torres et al. (2019) make use of a spatial scale threshold above which BMs dominate and below which ITs dominate. This spectral technique is less effective when the common spatial scale interval extends too much, which is particularly the case in winter due to the emergence of small vortices (<50 km) associated with mixed layer instabilities (Ajayi et al., 2020). Gonzalez-Haro et al. (2019) and Ponte et al. (2017) explore multi-sensor approaches with altimetry and sea surface temperature observations, motivated by the fact that BMs and ITs have distinct signature on both fields.
In this paper, we propose to simultaneously and dynamically map and separate SSH components of BMs and ITs from 2D altimetric observations, distributed in time. The proposed approach, illustrated on Figure 1, implements an alternating minimization algorithm that iteratively calls two mapping techniques, each one specific to the estimation of a component (BMs or ITs  Dimet & Talagrand, 1986). The algorithm is called 4DVar-SW hereafter. The joint estimation algorithm is tested with observing system simulation experiments (OSSE). The true state of the ocean (illustrated on Figure 2) is provided by a primitive-equation regional  10.1029/2021MS002613 4 of 17 model, ROMS, in an idealized configuration that simulates the propagation of a mode-one ITs through a turbulent mesoscale zonal jet Ponte et al., 2017).
The papers is structured as follows: Sections 2 and 3 present the data assimilation algorithms reconstructing BMs and ITs, respectively; Section 4 formalizes and describes the general joint estimation approach that implements both assimilation algorithms; The experimental set-up is presented in Section 5 and the results are discussed in Section 6; Conclusions and perspectives are drawn in Section 7.

Mapping the Balanced Motions: BFN-QG Algorithm
The mapping technique for BMs is presented in detail in Le Guillou et al. (2021). It is based on a QG model and a data assimilation method called Back and Forth Nudging (Auroux & Blum, 2008), and referred to as BFN-QG. The dynamical model is a 1.5-layer QG model (Ubelmann et al., 2015). Forcing and mixing terms are omitted, resulting in the simplified QG equations: where J is the Jacobian operator and the streamfunction ψ is proportional to SSH η: with g the gravity constant, and f the Coriolis parameter. The potential vorticity (PV) is linked to the streamfunction by the elliptical equation: where ∇ 2 is the Laplace operator and L R is the first baroclinic Rossby radius of deformation. Equation 1 needs boundary conditions (in practice, a prescribed SSH field is set at the border pixels; see Section 5 for more details). Propagating SSH from time t i to time t i+1 is a succession of elementary steps: (a) compute q i from η i with Equations 2 and 3, (b) propagate q i with Equation 1 to obtain q i+1 , (c) invert Equation 3 to retrieve ψ i+1 , then η i+1 .
The BFN technique is based on the Nudging method (Anthes, 1974), which consists of adding an extra term to the model prognostic equation (Equation 1) to pull the model variable towards the observations. This term is proportional to the difference between the observations and the model variable. Here, the model variable is PV q and the observations are 2D SSH η obs . Thanks to Equations 2 and 3, a corresponding PV observation q obs can be derived from η obs and used to nudge Equation 1 as: where K is a tunable coefficient controlling the strength of nudging. As explained in Le Guillou et al. (2021), the nudging term K(q obs − q) must exhibit smooth variations in time to avoid the emergence of numerical instabilities and enhance the constraint of the observations on the model dynamics. A Gaussian kernel is used for the nudging term: where K 0 and τ are the nudging parameters and t obs are the observation times. Note that, similarly, spatial smoothing needs to be implemented in the case of irregular space distribution of observations, as discussed in Le Guillou et al. (2021).
Nudging is also performed backward in time, by adding a feedback term with a sign opposite to that of the forward nudging. Note that the backward integration of the model is possible because the model Equation 1 is reversible in time.

of 17
The BFN is the combination of the forward nudging and the backward nudging over a sliding time window of length T BFN . Over a specific time window [t i , t i + T BFN ], the algorithm works as follows. The model state at time t i is prescribed from the previous time window. The forward nudging is first computed from t i to t i + T BFN . The final state at t i + T BFN is then used as initial condition to run the backward nudging back to t i . The result initializes the forward nudging, and so on until convergence. In a few iterations, the model converges towards a trajectory that both fits the observations and complies with the QG dynamics. The process can then be carried out in the next time window. Note that the time window T BFN must be long enough to include multiple observations and short enough to preserve dynamical consistency and ensure convergence within the window. The implementation details specific to the problem addressed in the paper are given in Section 5.

Mapping the Internal Tide Signal: 4Dvar-SW Algorithm
In this section, we aim at estimating the surface signature of ITs excited by one tidal component of frequency ω.

Dynamical Model
The dynamical model used to simulate the surface propagation of ITs motions is a linear SW model. This model represents the first baroclinic mode dynamics assumed to capture the largest part of the ITs signal (Ray & Zaron, 2016). The equations are (Gill, 1982): where (u, v) are the velocity components, η is the SSH, f is the Coriolis frequency, and H e the equivalent depth that determines the baroclinic deformation radius (through the relation = √ ∕ ). Equation 6 are discretized on a C-grid and a leap-frog time-stepping scheme (Sadourny, 1975). The domain is a square, land-free domain with a flat bottom and open boundaries. ITs enter through the boundaries and are not generated within the domain.
After entering the domain, ITs can be made incoherent by time variations of the equivalent depth, since H e modulates the phase speed of the waves, given by = √ 2 2 + ( ) 2 where k is the wavenumber. Space variations of H e (along with the spatial variation of f) changes the wave propagation direction and its wavelength (Rainville & Pinkel, 2006).
Open Boundary Conditions (OBCs) are designed to introduce and let out waves. They are of the radiative type (Blayo & Debreu, 2005). The amplitudes and directions of incoming waves are prescribed with a Flather condition (Flather, 1987) as: where v n is the normal component of the wave velocity at the boundary. Parameters and η ext must be prescribed. The signs in the above equation vary with the boundary (− for southern and western boundaries, + for the others).
Incoming waves are superpositions of progressive monochromatic waves with frequency ω and absolute wavenumber k. The frequency is a tidal frequency (details are given in Section 5) and k is deduced from ω using the dispersion relation: The superposing, monochromatic waves differ by their amplitudes and their propagation directions. Denoting and the amplitudes of velocity and SSH of the wave entering the domain with an angle θ with the inward 6 of 17 boundary normal direction, and denoting κ θ the unit vector of the propagation direction, the expressions of parameters and η ext at boundary grid point r = (x, y) and time t are: where i is the imaginary unit and is the complex notation.
Parameters and are linked together by the polarisation relations (obtained from Equation 6). Thus, prescribing boundary conditions with Equation 7 implies setting only at each boundary and for each angle θ. A limited number of angles are predefined. The prescription of the boundary conditions involves 2N θ N b N t values, with N θ the number of prescribed wave directions, N b the number of grid points at the boundary, and N t the number of model time steps. In the following, we denote all these external values to be prescribed by η ext .

4DVar Cost Function
The forward problem described in Section 3.1 is here formulated as: where = ( , ) , , are the input parameters, M is the non linear operator solving the prognostic Equation 6 associated with the OBCs (Equation 7) over a fixed time interval [t 0 , t 1 ]. In Equation 10, we neglect the initial condition in the function arguments. Given the high speed of the waves crossing the domain and the absence of balanced motions, the initial condition is completely forgotten after 10 days. As stated in Section 5, the 4DVar window is 4-month long, and the first (and last) 10 days are ignored for the performance diagnostics. The sea state essentially depends on boundary conditions η ext and model parameter H e . So, the initial condition is unimportant. In the following, we start experiments with a flat and motionless body of water (all variables to zero).
The associated inverse problem consists in seeking the set of parameters opt = ( opt , opt ) that produces a model state trajectory that best fits the observations η obs in the time window [t 0 , t 1 ]. We formulate this problem as the minimization of the 4DVar cost function: is the background control vector, H is the linear operator that projects the model state in the observation space, B and R denote the background and observation error covariance matrices, respectively. In this work, these matrices are set diagonal. This aspect is justified by the expression of the control variables in a reduced-order basis, making them independent from each other (see Section 3.3).
The minimum of J is found using a descent method that requires the calculation of the gradient ∇J with respect to ϕ: This gradient is computed with an adjoint method that makes use of the adjoint model M*. Since the model M is non linear, an incremental 4Dvar algorithm (Courtier et al., 1994) is implemented to find the optimal set of parameters ϕ opt that cancels ∇J.

Reduced-Order 4Dvar
Order reduction is a standard practice in geophysical data assimilation or inversion, to overcome the issues of ill-posedness and numerical complexity (e.g., Robert et al., 2005, for 4Dvar). Control variables are projected into a low-dimensional basis of vectors (EOFs -Empirical Orthogonal Functions-are a usual choice, relevant to many inverse problems) and the cost function minimization is performed in this low-dimensional space to make it possible or faster.

of 17
In this work, we propose to express the control vector (i.e., the model parameters ϕ) in a reduced basis G made of space-time Gaussian functions. This choice has been motivated by the fact that we know H e and η ext should vary "weakly" in space and time. Depth H e inherits its time variation scales from the balanced flow, and the tidal forcing is quasi-stationary in comparison with the generated internal tides. Mathematically, we call G (k,i) one basis function centered at the point (t k , r i ) (t k is a timestamp and r i is a spatial coordinate), expressed as: where τ (respectively σ) is the time (respectively space) characteristic scale. For H e , r and r i both refer to 2D coordinates (x, y) while for η ext they refer to a 1D coordinate (along the domain boundary). Thus, G (k,i) is a 3D function for H e and a 2D function for η ext . As an illustration, we represent in Figure 3 the space-time distribution of the basis functions for η ext at the western boundary and the wavy SSH produced by one element of the basis. Two successive basis functions are separated by t k+1 − t k = τ in time and ‖r i+1 − r i ‖ = σ in space.
In this basis, ϕ is written as: where w = (w ik ) are the coordinates of ϕ in the basis G, and the new control variables of the 4DVar. The cost function J and its gradient ∇J are then obtained by replacing H by H • G in Equations 11 and 12.
For the order reduction to be efficient, the characteristic scales τ and σ of the Gaussian basis functions have to be carefully chosen. They should be short enough to explain a maximum of the variability of H e and η ext and long enough to reduce significantly the dimension of the control space. Besides, the use of the order reduction (in addition to the motionless initial state) prevents non tidal dynamics (such as geostrophic motions) from appearing. This also impacts the choice of the characteristic scales: τ (respectively σ) must be larger than the tidal time period (respectively wavelength). Table 2 gives the retained (constant) values of τ and σ for the experiments presented in Section 5, considering the previous comments. With these values, the dimension of the control space is reduced by a factor of 10 4 . In addition to reducing the size of the control vector, the reduced order basis justifies the choice of a diagonal matrix for B, assuming that the basis element are independent from each other.

Joint Estimation Approach
As stated in the introduction, our goal is to map and separate BMs and ITs based on wide-swath altimetry observations that contains both components, and our approach is based on data assimilation. Let us represent by η BM the SSH field related to BMs, η IT the SSH field related to ITs, and η obs the observation of η BM + η IT . The reconstruction problem is formulated as finding: arg min where the norm typically considered here is the L 2 norm based on the space-time integral. The resolution of the minimization problem is addressed using a classical alternating minimization procedure (Tseng, 1990). The use of an alternating minimization algorithm is based on the strong assumption of the decoupling between the dynamics driving BMs and ITs. Iterations are performed, alternating the two minimization subproblems: In these equations, k is an iteration index. ̂ BM k and ̂ IT k are BMs and ITs estimations, respectively, at iteration k. Both minimization subproblems in Equations 16a and 16b can be solved with methods standard in geophysical data assimilation (4DVar, 3DVar, Kalman filters, etc). In what follows, Equation 16a is solved with BFN-QG (see Section 2) assimilating the corrected observation obs −̂ IT k-1 . Equation 16b is solved with 4DVar-SW (see Section 3) assimilating the corrected observation obs −̂ BM k . The joint estimation pseudo-code is presented in Algorithm 1 and illustrated in Figure 1.
With this setup, BMs and ITs dynamics are uncoupled (which is essential for the alternating assimilation to work): they are governed by distinct prognostic equations and no link is made between them (in particular, H e and L R are treated as independent quantities). In reality, BMs and ITs interact together. Here, we only simulate the effect of BMs on ITs propagation through the control of the space-time variation of H e .

Experiment Setup
The joint estimation algorithm is tested with OSSEs: a simulation provides reference surface data containing both BMs and ITs motions; SSH observations are extracted from this simulation and are assimilated; the estimated fields are then compared to the reference fields. Algorithm 1. The joint estimation algorithm. Variables ̂ BM and ̂ IT are BMs and ITs estimations, respectively. C is a convergence criterion. In this study, we stop the algorithm when the space-time trajectories of BMs and ITs do not differ more than 1% of those obtained at the previous iteration.

Reference Simulation
The test flow consists of a turbulent zonal jet and a mode-1 internal tide propagating signal described in Dunphy et al. (2017) and Ponte et al. (2017). This simulation is carried out with the Coastal and Regional Ocean COmmunity (CROCO, https://www.croco-ocean.org) model. The horizontal resolution is 4 km and the domain size is 1,024 km × 2,880 km. The mode-one internal tide is generated by a zonally uniform tide excitation in the southern part of the domain. The time period of the tide is T = 12 hr. The plane wave produced is propagating northward and is encountering the meridional jet located at the center of the domain (Figure 2).
The reference SSH, called η truth , is decomposed into BMs and ITs . The reference BMs field BM truth is obtained by low pass filtering η truth : . The reference ITs field IT truth , is obtained with an harmonic fitting:  (Zaron & Egbert, 2014) for the ITs. The ITs wavelength is about 150 km. South of the jet, the wave is rectilinear and aligned with the southern boundary, which reflects its stationarity. North of the jet, the wave is dispersed due to its interaction with turbulent BMs, producing nonstationarity in the wavefield . The (weak) residual of this decomposition (| truth − BM truth − IT truth | < 1 cm, not shown) corresponds mostly to the mode 2 ITs generated by the interaction of the wavefield with the turbulent jet.
In the following experiments, the horizontal resolution is degraded to 16 km (by local averaging) to reduce computational complexity. Observations are extracted from the degraded fields. Both models (QG/SW) and diagnostics are computed on this coarse grid.

Observational Scenarios
Observations are drawn from the true (degraded) SSH (η truth ), containing both BMs and ITs components. They are free of error, complete in space, and regularly distributed in time. The observation frequency is 75 hr in the reference experiments discussed in Section 6.1.
In order to help evaluation of the mapping performance of the joint estimation experiments (hereafter denoted as joint), two additional scenarios are considered. The first one, denoted idealized, refers to the case where both components can be observed separately, and each observation is processed with the appropriate assimilation system. This means that BFN-QG (4Dvar-SW) is fed with observations extracted from BM truth ( IT truth , respectively) fields. This ideal scenario provides upper bounds for the mapping performance of the joint experiments. The second scenario, hereafter denoted as independent, implements each assimilation system independently (without iteration) with observations of the full signal η obs . This scenario provides a baseline to be outperformed by the joint experiments.
This study also investigates the impact of the observational temporal sampling on the mapping performances. For that purpose, a sensitivity experiment is performed and discussed in Section 6.2.

Configuration of Assimilation
Assimilation is performed in an inner domain of the reference simulation, illustrated by the black squares in Figure 2, to avoid the effects of boundary conditions and forcing in the reference simulation. The experiments are run over 4 months.
The BFN-QG parameters are listed in Table 1. The Rossby radius of deformation L R is chosen from a climatology for western boundary currents. Note that L R is fixed, while H e is a space-time varying control variable of the 4DVar-SW system. The boundary conditions for the QG model were given by the mean SSH field, which is computed by time averaging the reference SSH η truth over the entire simulation time period. For operational implementation, DUACS L4 product (or equivalent) could be used, as in Le Guillou et al. (2021). The length of the BFN time window T BFN is fixed to 15 days after a few trial and error experiments. After BFN convergence (reached when the forward and backward trajectories do not differ from the previous ones by more than 1%), the 7 days in the middle of the interval (days 5-11) are saved and the BFN windows is shifted by 3.5 days for the next iterations. Results from days shared by two consecutive windows are averaged. This guarantees smoothness in the time sequence of SSH.
The 4Dvar-SW parameters are listed in Table 2. For the joint scenario, the background values b and b are updated at each iteration by the optimal values found at the previous iteration. We consider two values for the observation error covariance matrix R: 0.01 m for the joint and idealized scenarios, 0.1 m for the independent scenario. The value for the independent scenario is higher to reflect representativeness errors and allow convergence of the 4Dvar-SW even with the presence of BMs in the observations. We set one order of magnitude below the amplitude of the target signal because of the superposition of the Gaussian basis functions in space (as illustrated in panel b of Figure 3).

Diagnostics
For measuring the performance of the reconstruction, we use the time varying Root Mean Square Error (RMSE) and its associated score S, defined as follows: where N is the number of pixel in the study domain, ̄ is the time-average of the variable X and RMS is the root mean square function. A score of 1 indicates a perfect reconstruction in terms of RMSE, while a score of 0 indicates that the time-averaged RMSE is as large as the RMS of the reference field. Note that for BMs, the score is defined relative to the time anomaly of the reference field. Figure 4 illustrates the performances of the joint estimation algorithm over the iterations for a temporal sampling of 75 hr (6.25 ITs periods). The joint estimation algorithm converges after 10 iterations. The scores obtained with the idealized and independent scenarios are also shown. As the joint estimation is started by a BMs reconstruction, the first iteration of the joint algorithm corresponds to the independent scenario for BMs. Without the alternating minimization, the estimation scores for both BMs and ITs fall between 60% (worst case, independent scenario) and 80% (best case, idealized scenario). The 80% bound is likely due to the observation time sampling and to physical processes unresolved by the QG and SW models. These processes include the presence of higher baroclinic modes and ageostrophic dynamics for the QG model, higher tidal harmonics and neglected terms in the SW model. The lower 60% bound reflects the presence of the component that the system cannot (or hardly) reconstruct (ITs for the QG model, BMs for the SW model). This component appears as noise in the observations. The joint estimation algorithm is able to disentangle BMs and ITs components of SSH. Throughout iterations, both components are progressively separated. This is particularly visible in the BM snapshots, where ITs signature is clearly visible at the first iteration, and gradually filtered out throughout iterations. We attribute the strength of the joint estimation to the fact that each data assimilation sub-system (BFN-QG and 4DVar-SW) filters out the component it is not supposed to reconstruct. This is reflected in the relatively high (closer to 1 than 0) scores reached in the independent scenarios.

Reconstruction Performances
The 4Dvar-SW technique can reconstruct a large part of the nonstationary part of ITs. Figure 5 shows snapshots of the stationary and nonstationary parts of the reference and estimated ITs fields. Stationary part has been extracted with an harmonic fitting of the signal over all the experimental window. Nonstationary part is obtained by subtracting the stationary part from the full ITs field. Eight seven percent of the stationary and 44% of the nonstationary ITs variances are recovered on average. For the nonstationary ITs fields, 4Dvar-SW captures 54% of the variance in the northern part of the domain, but only 28% in the southern part. This poor performance is attributed to the weakness of the nonstationary signal in the southern part ( Figure 5, right column). Besides, some physical processes, such as ITs reflection on the turbulent jet, refraction by relative vorticity gradients and vertical shear Olbers, 1981) may not be captured by 4Dvar-SW.
The spatial patterns of SSH BM estimated by BFN-QG and the optimal equivalent depth opt obtained with 4Dvar-SW share similarities together and with the reference equivalent depth, as illustrated on Figure 6. A natural improvement of the method would be to use BM to estimate a background value b in the 4Dvar-SW to improve the ITs reconstruction. This would require to exploit the fact that H e reflects the "large scale" baroclinic radius (L R ) modified by the mesoscale eddies through modification of the stratification profile N(z) Kelly & Lermusiaux, 2016). However, to do so, one would need a climatology of the stratification and, more critically, a knowledge of the vertical structure of the balanced flow, in order to convert SSH anomalies to H e . Reciprocally, 4DVar-SW could help refine the estimation of H e and stratification, thus providing some information on the vertical structure of the balanced currents and on the ocean heat content (Zhao, 2016).

Sensitivity to the Observation Temporal Sampling
In this section we investigate the impact of the observation temporal sampling on the estimation of ITs and BMs. Thirteen time intervals Δt obs are tested, from 69 to 75 hr by increments of 30 min. The ITs period being T = 12 hr, the 6 hr-interval length corresponds to T/2 and the central 72 hr step corresponds to a multiple of T. For each Δt obs , we run an idealized, an independent and an joint experiments. The results are shown on Figure 7.
The ITs estimation is strongly degraded near Δt obs = 72 hr, which affects the reconstruction of BMs in the joint estimation. At Δt obs = 72 hr, ITs estimation scores are even negative in the joint and independent scenarios. This particular time interval presents a singularity and leads to an aliasing issue, since it observes only one phase of the wavefield and does not detect the wave propagation. The same occurs with multiples of T/2 (e.g., 66 and 78 hr, not shown here), when only two opposite phases of the wavefield are observed. Since the ITs component is not filtered out by 4DVar-SW, the estimation of BMs by BFN-QG in the joint scenario is similar to the independent one.
Using the joint estimation algorithm mitigates the poor ITs reconstruction with time intervals near 72 hr. The ITs estimation scores with Δt obs = 71 hr or 73 hr are negative in the independent experiments, but exceed 0.6 in the joint experiments. This experiment suggests that specific temporal samplings of wide-swath altimeters (e.g., that

Conclusions
In this study, we have presented an alternating minimization algorithm to simultaneously and dynamically map and separate ITs and BMs signals from 2D altimetric observations. Two data assimilation techniques are used in an iterative process to estimate both signals: BMs are reconstructed by the BFN-QG algorithm while ITs are reconstructed by the 4Dvar-SW algorithm. At any given iteration, the assimilated observations used to map one signal are made of the difference between the full observation (containing both signals) and the previous estimation of the other signal. The use of the 4Dvar-SW algorithm for ITs reconstruction allows to map (i.e., estimate in magnitude and phase) both stationary and nonstationary ITs signals which is not obvious with other conventional separation methods such as harmonic analyses. The joint estimation of BMs and ITs has been tested in an idealized context, with a true state of the ocean and artificial observations extracted from a numerical simulation that propagates a mode-one internal tide through a turbulent quasigeostrophic jet. This numerical experiment aims at illustrating the feasibility of this approach in a simplified context and should constitute a step closer to the processing of wide-swath altimetry data.
The results show that the joint estimation allows to simultaneously separate and estimate 78% and 75% of the BMs and ITs variance, respectively. The remaining variance is interpreted as being due to non resolved dynamics by the QG and SW models.
A sensitivity study of the observational temporal sampling has also been performed and shows that an unfortunate temporal sampling of the observations can be detrimental to the joint estimation performances. In particular, the performances of the ITs reconstruction are considerably reduced when the time step between two consecutive observations is a multiple of half the ITs time period T/2. Fortunately, it turns out that this situation is singular for non sun-synchronous orbits, meaning that a slight shift from a multiple of T/2 in the observation time sampling ensures good performances again.
Although the results of the observational temporal sampling sensitivity study are encouraging, future work should investigate the method applied to a realistic observational sampling, such as SWOT. One of the major next steps will be to assess the potential of the method with spatially sparse observations. An interesting avenue to improve the method in that context could be to use additional information, such as conventional nadir data, to ensure the joint estimation method convergence hence improving the reconstructions.
Another crucial step towards operational SSH mapping will be to experiment on a more realistic ocean simulation, with multiple tidal components and harmonics in constant interaction. To allow the 4Dvar-SW to efficiently extract ITs signal from realistic SSH data, some technical improvement may be needed. Source terms, that generate ITs signal inside the study domain, may have to be implemented and controlled by the assimilation to map ITs signals around generation sites (located at ocean ridge and sea mounts for instance). One can also think of improving the resolved dynamics by including higher baroclinic modes. Finally, it may be needed to add interaction terms between BMs and ITs in the shallow water equations (as described in Kelly & Lermusiaux, 2016) to correctly estimate complex ITs field. For that last remark, one first strategy would be to estimate a background term of H e from the reconstruction of BMs.