When evaluated against the 1/30°-resolution, submesoscale-resolving OFES model outputs, the previously published “interior + surface quasigeostrophic” method (from the 2013 study by Wang et al., denoted W13) for reconstructing the ocean interior from sea surface information is found to perform improperly in depicting smaller-scale oceanic motions (associated with horizontal scales smaller than about 150 km). This could be attributed to the fact that the W13 method uses only the barotropic and first baroclinic modes for the downward projection of sea surface height (SSH), while SSH at smaller scales significantly reflects other higher-order modes. To overcome this limitation of W13, an extended method (denoted L19) is proposed by employing a scale-dependent vertical projection of SSH. Specifically, the L19 method makes the projection via two gravest modes as proposed in the W13 method only for larger-scale (>150 km) signals, but for smaller scales (≤150 km) it exploits the framework of the “effective” surface quasigeostrophic (eSQG) method. Evaluation of the W13, eSQG, and L19 methods shows that the proposed L19 method can achieve the most satisfactory subsurface reconstruction in terms of both the flow and density fields in the upper 1000 m. Our encouraging results highlight the potential applicability of L19 method to the high-resolution SSH data from the upcoming Surface Water and Ocean Topography (SWOT) satellite mission.
Over the past decades, advancements in satellite technology have revolutionized oceanography by providing global-scale measurements of sea surface variables. The future Surface Water and Ocean Topography (SWOT) satellite mission aims at capturing, on a global average, sea surface height (SSH) signals with wavelengths ≥ 15 km (Durand et al. 2010; Fu and Ubelmann 2014; Wang et al. 2019) and is expected to further improve our understanding of the oceanic mesoscale (from 300 down to 50 km) and submesoscale (<50 km; Su et al. 2018; Torres et al. 2018) processes. However, satellites only capture the sea surface information. The scarcity of subsurface observations makes the method for reconstructing the ocean interior from satellite measurements a key issue in oceanography, especially the method geared toward the upcoming high-resolution SWOT SSH.
While statistical methods relying on sufficient historical statistics can achieve quite satisfactory reconstructions (see Chu et al. 2014; Fresnay et al. 2018, and references therein), dynamical methods (e.g., Cooper and Haines 1996; Zhang et al. 2014) can obviate the availability of sufficient historical data and provide dynamically consistent fields. In recent years, the surface quasigeostrophic (SQG) theory, which can reasonably describe the dynamics induced by the sea surface density (SSD) acting as a surface potential vorticity (PV) sheet (Bretherton 1966), has been explored for the subsurface reconstruction by many investigators (see Lapeyre 2017, and references therein). Lapeyre and Klein (2006) adjusted the SQG theory and proposed an “effective” SQG (eSQG) approach to reconstruct the upper-ocean structures from either the SSH or SSD. This approach takes into account the nonnegligible contribution from PV in the ocean interior (not involved in the SQG formalism), and has been successfully employed in several studies (e.g., LaCasce and Mahadevan 2006; Isern-Fontanet et al. 2006, 2008; Klein et al. 2009; LaCasce 2012; Ponte and Klein 2013; Smith and Vanneste 2013). Utilizing SWOT-like SSH data as input of the eSQG approach, Qiu et al. (2016) obtained satisfactory reconstruction of three-dimensional circulation (relative vorticity and vertical velocity) fields in the upper 1000 m, which highlights the potential applicability of this approach to the future SWOT measurement. Although the eSQG method has been demonstrated to be promising, its underlying assumption of good correlation between SSD and interior PV (in which case SSD is in phase with SSH) may lead to degraded reconstructions when there exists a noticeable phase shift between SSD and SSH (Wang et al. 2013; González-Haro and Isern-Fontanet 2014; Liu et al. 2014; Isern-Fontanet et al. 2014, 2017).
Wang et al. (2013) proposed a more elaborate method termed “interior + surface QG” (hereinafter W13) to reconstruct the subsurface structures from both SSD and SSH. This method first determines the dynamical solution associated with surface PV (i.e., SSD) via the SQG theory; it then determines the solution induced by interior PV, through two gravest vertical normal modes, whose amplitudes are deduced from the residual SSH (the difference between the total SSH and the SQG contribution). In contrast to the eSQG approach, the W13 method utilizes both SSD and SSH to simultaneously constrain the surface and interior PV, and does not require good correlation between SSD and SSH. The W13 method has been validated using fields from numerical simulations (Wang et al. 2013; Liu et al. 2014; LaCasce and Wang 2015), in situ and satellite observations (Liu et al. 2017), all of which yielded promising results.
A satellite-derived two-dimensional SSH map, whose spatial resolution is currently confined to wavelengths longer than about 150 km (Chavanne and Klein 2010; Chelton et al. 2011; Dibarboure et al. 2011; Gaultier et al. 2016; Wang et al. 2018; Wang and Fu 2019), mainly reflects the barotropic (BT) and first baroclinic (BC1) modes in the oceans (Stammer 1997; Wunsch 1997, 2015; Smith and Vallis 2001). Based on this argument, the W13 method projects the residual SSH downward through BT and BC1 to determine the interior PV contribution. However, as indicated by the results of Badin (2014) and LaCasce (2017), including only BT and BC1 is insufficient for capturing the smaller-scale SSH signals (unresolved by the present satellite product) that are significantly composed of higher-order modes in addition to BT and BC1 (Lapeyre 2009; Wunsch 2015). Therefore, the two-gravest-mode scheme imposes a limitation on the W13 method: the method cannot be properly applied to the high-resolution surface information (such as SWOT SSH). Such a limitation has not been thoroughly investigated by previous studies.
There are two goals for this article. One is to investigate the limitation of the W13 method in reconstructing the ocean interior from high-resolution sea surface data. The second is to extend the W13 method to overcome the limitation, by improving the vertical projection scheme for the residual SSH. The rest of this paper is structured as follows. The 1/30°-resolution Ocean General Circulation Model for the Earth Simulator (OFES) simulation that serves as the testbed of our explorations is described in the next section. Section 3 briefly reviews the W13 method and then expounds its limitation. In section 4 we extend the W13 method by taking the advantage of the eSQG approach, and evaluate the performance of the extended method (hereinafter L19). Section 5 summarizes the results from the present article and provides perspectives for the future work.
a. The 1/30°-resolution OFES simulation
In the present study, we use the OFES simulation (Masumoto et al. 2004; Komori et al. 2005) of the North Pacific Ocean (20°S–66°N and 100°E–70°W) at high resolution in both the horizontal (1/30°) and vertical (100 levels) directions as “true” fields to evaluate the reconstruction methods. This 1/30° simulation uses outputs from the global 1/10° OFES simulation (Nonaka et al. 2016) on 1 January 2000 as the initial condition and is forced by the long-term surface wind stresses and heat flux data of the 6-hourly Japanese 25-year Reanalysis (JRA-25) product (Onogi et al. 2007). For horizontal mixing of momentum and tracers, the model uses a biharmonic operator to reduce numerical noise. For vertical mixing, the scheme developed by Noh and Kim (1999) is employed. Details of the model can be found in Sasaki et al. (2014, 2017). Notice that, as pointed out by Qiu et al. (2016), the 1/30° simulation can capture the variability of length scale similar to that expected from the SWOT satellite.
The 1/30° OFES simulation was integrated for 3 years (Sasaki and Klein 2012), and the daily averaged data of temperature, salinity, horizontal velocity, and SSH fields in the last two years (i.e., 2001 and 2002) are utilized in our following analyses. Notice that the daily average filters out a large part of near-inertial motions, and only the balanced part of the flow is retained (Sasaki et al. 2017). We conduct the upper-ocean reconstruction experiments using the surface information from the OFES simulation as input, in the subdomain of the Kuroshio Extension with a size of 10° × 10° (dashed box in Fig. 1). Reconstructed flow and density fields are assessed against the simulation itself in a smaller 6° × 6° region (solid box in Fig. 1) to avoid the edge effect (Qiu et al. 2016).
b. Data preprocessing
The input of reconstruction methods includes three variables: sea surface density anomaly field (SSDA), sea surface height anomaly field (SSHA), and vertical stratification profile N2(z).
Based on the UNESCO1983 equation of state (Fofonoff and Millard 1983), we first derive the potential density fields at different levels using temperature and salinity data from the OFES simulation. Then the potential density fields are used to calculate an area-averaged N2(z) utilizing the mean pressure between two vertical grids as the reference pressure:
Here ρυ and are respectively the volumetric and horizontal means of potential density in the target region (dashed box in Fig. 1). Such N2(z) could not be directly used in the reconstruction because it approaches zero in the shallow layers, which can lead to unrealistic overshoots (Wang et al. 2013). Therefore, we make adjustments to the N2(z): replacing the surface value by a mixed layer–averaged one; then a linear interpolation between this new surface value and the value at the base of mixed layer is employed as the adjusted stratification profile in the mixed layer (Wang et al. 2013; Liu et al. 2014, 2017). A constant stratification is also used in this study (see section 4). According to Isern-Fontanet et al. (2008, 2017), we set to be an averaged value of the N2(z) (derived directly from OFES instead of the vertically adjusted profile described above) in the top 1000 m.
To get the SSD and SSH anomalies, we remove the respective large-scale background information from the SSD and SSH fields. Background information is derived through a least squares fit of a field to the quadratic surface S(x, y) (Wang et al. 2013; Liu et al. 2014),
Here x and y represent the zonal and meridional coordinates, respectively. We then apply the two-dimensional trapezoid window to the anomaly fields in a 1° band along the box edges (for the purpose of Fourier transform), the same way as Qiu et al. (2016) did.
As an example, Fig. 2 presents the SSDA, SSHA, and stratification profiles on 31 March 2001 (corresponding to the area limited by the dashed box in Fig. 1). Note that, as a consequence of the SQG theory (Klein et al. 2008; Isern-Fontanet et al. 2008; Sasaki and Klein 2012; Lapeyre 2017), the SSDA field exhibits smaller scales than SSHA. In Fig. 2c, dotted line represents the N2(z) directly derived from the OFES simulation; the black solid line denotes the adjusted N2(z) (with an area-averaged mixed layer depth being 100 m); and gray line shows the constant .
3. Application of the W13 method to the OFES output
In this section, we use the OFES fields on 31 March 2001 (results are typical for all days throughout the 2-yr OFES simulation period) to evaluate the W13 method, which is briefly described in the appendix [a more detailed description can be found in Wang et al. (2013)]. The W13 method first determines the SQG component of the streamfunction Ψsur(z) using SSDA (Fig. 2a) and stratification (black solid line in Fig. 2c), and obtains the residual SSHA (the difference between the total SSHA and the SQG-SSHA, see Figs. 2b and 3a,b). To determine the interior component of the streamfunction Ψint(z), the method then projects downward the residual SSHA [i.e., Ψint(z = 0)] using BT and BC1 with the no-slip boundary condition, that is, zero velocity, at the bottom.
We focus on three variables deduced from the reconstructed streamfunction Ψ: horizontal velocity anomaly V = z × ∇Ψ, relative vertical vorticity anomaly ζ = ∇2Ψ, and density anomaly ρ = −(ρ0f0/g)(∂Ψ/∂z).
a. Limitation in flow field reconstruction
Figure 4a shows the horizontal velocity anomaly vectors from the W13 method (blue arrows) and OFES simulation (red arrows) at the 500-m depth. Although the W13 solution is satisfactory with correlation between the reconstruction and OFES data being 0.89 (0.9) for the zonal (meridional) velocity, a close scrutiny indicates that the method cannot properly depict the smaller-scale structures (such as those around 35°N, 144.5°E). Discrepancies between the two velocity fields become more notable in the deeper layers (not shown).
A more stringent test comes with the comparison of vorticity ζ, as this highlights small horizontal scales in the flow field (Lapeyre 2009; Wang et al. 2013). Plane views of 500-m ζ fields clearly display the W13 limitation: the W13-reconstructed ζW13 (Fig. 5b) almost fails to capture the OFES simulation ζOFES (Fig. 5a), with the pattern correlation between them being only 0.42. A closer look at Figs. 5a and 5b reveals that while similarities at larger scales are obvious, it is the spurious smaller-scale structures in ζW13 (Fig. 5b) that lead to the low pattern correlation. This can be more clearly illustrated in the spectral space. At wavelengths shorter than about 150 km the ζW13 spectral density (black line in Fig. 6a) unreasonably exceeds the OFES value (gray line in Fig. 6a); meanwhile, the spectral correlations between ζW13 and ζOFES are low (black line in Fig. 6b), especially at scales smaller than 100 km. It is important to note that the W13 method can still achieve a satisfactory reconstruction at scales larger than about 150 km.
The abovementioned spurious smaller scales (<~150 km) exist in the ζW13 field at all depths (not shown) below the mixed layer (100-m depth), which can also be demonstrated in the spectral space: the ζW13 spectrum (Fig. 7b) is much larger than the observation (Fig. 7a) and the spectral correlations are weak (Fig. 8a). These spurious signals notably degrade the W13 performance in physical space: pattern correlations between ζW13 and ζOFES, which have a decreasing tendency with depth in the 100–1000-m layers, fall below 0.5 beneath the 400-m depth (black line in Fig. 9a). Notice that the W13 limitation hardly arises in the mixed layer where the W13 reconstruction resembles the true field quite well with pattern correlations exceeding 0.8 (black line in Fig. 9a).
Decomposition of the 500-m ζW13 (Fig. 5b) into the SQG (ζsur, Fig. 5d) and interior (, Fig. 5e) components shows that the spurious smaller scales can be attributed to the interior portion determined through the downward projection of the residual SSHA (Fig. 3b). In other words, these smaller scales originate from the residual SSHA [i.e., Ψint(z = 0)], as also indicated by the high similarity (with pattern correlation reaching 0.99) between the residual-SSHA-derived ζ (Fig. 3c) and the 500-m (Fig. 5e). In the following section 3b, we will specifically analyze the reason for the W13 limitation arising from the improper SSHA downward projection.
b. Deficiency of only BT and BC1
Indeed, as the submesoscale features (<50 km) are mostly related to mixed layer instabilities (Sasaki et al. 2014), such features cannot be explained by the traditional Phillips (1954) model of baroclinic instability that influences larger-scale dynamics and that rationalizes the BT/BC1 representation (Klein et al. 2008). Meanwhile at the short mesoscale range (50–150 km; Li et al. 2019), some surface-intensified structures also do not pertain to the middepth-intensified Phillips-type instability (Smith 2007; Tulloch et al. 2011; Rocha et al. 2013). To expound the W13 limitation, we remove the SQG contribution from OFES data to get the residual velocity (anomaly) fields, based on which the residual kinetic energy (KE) is calculated. Then we examine the partition of residual KE into BT/BC modes. A caveat (Lapeyre 2009) is that the sum of the residual KE and SQG-KE differs from the total (OFES) KE, since the SQG mode and BT/BC modes are not orthogonal (Wang et al. 2013).
Figure 10a shows the surface KE spectra of the different modes (Fig. 10b). Most notably, the BT importance (black solid line) decays fast as the scale decreases, while BC2 and BC3 (gray lines) become as important as BT for wavelengths close to 150 km and contribute more than BT at wavelengths shorter than about 100 km (where other higher-order modes are also significant, not plotted for clarity). Although BC1 (black dashed line) indeed dominates for almost all wavelengths shorter than 600 km, Fig. 10a indicates that the W13 method employing only BT and BC1 cannot properly project downward the smaller-scale (<~150 km) residual SSHA [or, Ψint(z = 0)] notably reflecting higher modes. Figures 10c–f (red line) further examine the root-mean-square (rms) of residual KE at different depths in different wave bands. For smaller horizontal scales (<~150 km, Figs. 10c–e), the BT partition (black solid line) accounts for a small fraction of the variance (11%, 16%, and 19% for the <50-, 50–100-, and 100–150-km wave bands, respectively); meanwhile, a larger fraction of the variance is in BC2/BC3 (gray lines) and higher modes (43%, 31%, and 23% for the <50-, 50–100-, and 100–150-km bands). Accordingly, these smaller scales (<~150 km) decay faster with depth than the larger ones (150–300 km), whose partition into BT (32%) and BC1 (57%) is prominent (Fig. 10f). These results suggest that the deficiency of only BT and BC1 in the W13 framework can overly project Ψint(z = 0) (i.e., residual SSHA), as the two gravest modes exert much greater influences to greater depths (especially BT that is depth invariant). As a result, smaller-scale features in Ψint(z = 0) do not attenuate fast enough and manifest falsely as spurious signals in the deeper layers. In the case of vorticity reconstruction (section 3a), the second-order horizontal derivatives highlight smaller horizontal scales and thus can make those spurious signals in Ψint(z) more prominent (cf. ζOFES and ζW13 spectra in Figs. 7a,b; see also the low spectral correlations in Fig. 8a).
It is important to notice that, at the 150–300-km mesoscale wavelengths (Fig. 10f), BT and BC1 projections (black lines) are quite dominant, and these two account for 89% of the variance. This rationalizes the two-gravest-mode scheme in the W13 method that still can perform satisfactorily at these scales (see Figs. 7a,b and 8a).
c. Density field reconstruction
The W13 method hardly manifests its limitation in the density reconstruction. As shown by Fig. 11 at the 500-m depth, the combination of the SQG (ρsur, Fig. 11d) and interior (, Fig. 11e) components leads to a successful W13 solution ρW13 (Fig. 11b) capturing the OFES simulation ρOFES (Fig. 11a) quite well, with the pattern correlation reaching as high as 0.95. Unlike the case of ζ reconstruction (Fig. 5e), here the interior component (Fig. 11e) is satisfactory. This is because the W13-reconstructed smaller-scale Ψint signals that overly project onto BT (and thus falsely penetrate into deeper layers) are eliminated in the density field, which is related to the vertical derivative of the streamfunction. As also suggested in the spectral space, at smaller scales (<~150 km), the ρW13 spectrum (Fig. 12b) is similar to that of the true field from OFES (Fig. 12a) in the upper 1000 m (cf. black and gray lines in Fig. 12d).
The W13 density reconstruction is successful throughout the upper 1000 m. Although the spectral correlations between ρW13 and ρOFES at shorter wavelengths (<150 km, especially <~100 km) are low in the deep layers (Fig. 13a; black line in Fig. 13c), they hardly affect the total field in the physical space (Qiu et al. 2016) where the pattern correlations exceed 0.9 (black line in Fig. 9c). Notice that at 0 m the pattern correlation presents a value of 1.0, since the W13 density solution at the surface is identical to the SSDA from OFES.
4. The L19 method and its evaluation
a. Proposition of the L19 method
The key point of the L19 method is that it introduces a cutoff scale LC and makes a scale-dependent, downward projection of SSHA to deduce Ψint. As suggested by Wortham and Wunsch (2014), the partition of SSHA into BT and BC structures as a function of horizontal scale could be the greatest (theoretical and observational) challenge in oceanography. According to discussions in the introduction and analyses in section 3b, we preliminarily set the cutoff scale LC to be 150 km in this study. For larger horizontal scales (>LC), at which the BT and BC1 contributions are dominant (see Fig. 10) and the W13 method can perform satisfactorily (see Figs. 7a,b and 8a), we retain the W13 two-gravest-mode scheme for the residual SSHA projection [i.e., Eqs. (A4)–(A6) in the appendix].
For smaller scales (≤LC) at which BC2 and other higher-order modes are more important (see Fig. 10), we choose to exploit the eSQG framework (see the appendix for details) that can work effectively (e.g., Klein et al. 2009; Ponte et al. 2013; Chavanne and Klein 2016). Specifically, we find that the difference between Eqs. (A9) and (A8) indicates structures of Ψint decaying away from the surface exponentially at a rate determined by their scale (the smaller scale decays faster than the larger one) and a constant stratification N0: exp[(N0/f0)κz]. Therefore, we employ this eSQG-derived vertical structure to deduce Ψint, whose amplitude is determined by matching at the surface to the residual SSHA [implicitly defined in Eq. (A5), i.e., ]:
Now the interior solution Ψint at all scales can be determined, and thus the total solution Ψ.
As Lapeyre and Klein (2006) have suggested, in the first 500 m, the surface density generally gives a larger contribution than the interior PV. Therefore Ψint only modifies the vertical variation of Ψ (Isern-Fontanet et al. 2008), which rationalizes the approximation of Ψint by the SQG-like solution in Eq. (3). It is important to notice that at smaller scales (≤LC), the SQG-like exponential structure exp[(N0/f0)κz] can reasonably capture the higher-order BC modes, as demonstrated by LaCasce (2012) in investigating how the SQG-like motion would appear in terms of the vertical normal modes. Therefore, the proposed L19 method can circumvent the W13 limitation arising from the deficiency of only BT and BC1. In fact, the SQG-like solution and traditional flat-bottom QG modes do not exclude each other (Rocha et al. 2013): at long wavelengths, the SQG solution quite resembles the BT/BC1 combination. For an ocean with rough topography, the BT mode vanishes (in reality, it is replaced by a bottom-intensified topographic wave mode; see LaCasce 2012), and the similarity between the SQG solution and BC1 is even more striking (de La Lama et al. 2016).
In what follows, we first demonstrate the L19 reconstruction on 31 March 2001, and then analyze the statistical results over the two simulation years (i.e., 2001 and 2002). Notice that for the horizontal scales smaller than LC, the L19 method uses a constant stratification profile (e.g., gray line in Fig. 2c).
b. Comparison with the W13 reconstruction
For flow fields at the 500-m depth, the L19 method (blue arrows in Fig. 4b) better represents the observed horizontal velocities (red arrows) than the W13 method (blue arrows in Fig. 4a), with pattern correlation for the L19 solution reaching 0.93 (0.95) for the zonal (meridional) velocity (see Table 1). Furthermore, the L19-reconstructed vorticity ζL19 (Fig. 5c) noticeably outperforms ζW13 (Fig. 5b): ζL19 satisfactorily captures ζOFES (Fig. 5a) with the pattern correlation reaching 0.75 (Table 1). The L19 method performs better because it can achieve a reasonable interior component (Fig. 5f), as those spurious smaller scales in (Fig. 5e) no longer emerge in .
The advantage of L19 over W13 can be more clearly illustrated in the spectral space (see Figs. 7 and 8 for the wavelengths shorter than 150 km). Below the mixed layer (100-m depth), the ζL19 spectrum (Fig. 7c) closely resembles the SQG-like structures of the OFES simulation (with the smaller horizontal scale decaying faster than the larger one in the vertical, see Fig. 7a) and is evidently better than the ζW13 counterpart (Fig. 7b), although the ζL19 spectral density is a little smaller than the actual value (cf. red and gray lines in Fig. 6a). Meanwhile, the spectral correlations for ζL19 (Fig. 8b) exceed those for ζW13 (Fig. 8a), noticeably at the short mesoscale 50–100-km band (also compare red and black lines in Fig. 6b). It is important to note that at submesoscales (<50 km), the L19 method (or the eSQG framework acting as a low-pass filter and damping smaller scales for a given depth) can only represent the tendency to diminish with depth but cannot fully reproduce these submesoscales (see the upper 100-m layer in Fig. 8b), as also suggested by Isern-Fontanet et al. (2008) and Qiu et al. (2016).
As seen in Fig. 9a showing pattern correlations between the reconstructed and OFES ζ fields (on 31 March 2001), the L19 method (red line) outperforms the W13 method (black line) below the mixed layer (100-m depth), with correlations for ζL19 falling in between 0.59 and 0.88 in the upper 1000 m. We further obtain the statistical results by averaging the daily pattern correlations over two years, which are presented in Fig. 9b. Statistics confirm that the L19 method can successfully reconstruct the ζ field in the 1000-m upper ocean: pattern correlations (red line) exceed 0.8 in the first 250 m, and fall in between 0.66 and 0.8 below the 250-m depth. In contrary, the W13 performance (black line) is poor below the 480-m depth (with pattern correlations being lower than 0.5), especially in the deeper layers. Notice that in the upper 50 m where the W13 limitation hardly arises, the two methods perform similarly well.
Now we compare the ρ reconstruction. Below the mixed layer, as analyzed in section 3c, the spurious smaller scales in W13-reconstructed Ψint (caused by the false partition into BT/BC1) have been greatly eliminated by the vertical derivative. Therefore, the advantage of L19 over W13 is inconspicuous. As shown at the 500-m depth, the L19 interior component (Fig. 11f) and total solution ρL19 (Fig. 11c) are almost the same as their W13 counterparts (Figs. 11e,b), with pattern correlation between ρL19 and ρOFES (Fig. 11a) still being 0.95 (Table 1). This is also indicated by Figs. 12c and 12b: ρL19 and ρW13 present quite similar spectral density below the 100-m depth (also see red and black lines in Fig. 12d). Indeed, below the base of mixed layer, pattern correlations (in the physical space, see Figs. 9c,d) between ρL19 and ρOFES (red lines) are almost identical to those for ρW13 (black lines), although the L19 method actually achieves better spectral correlations at smaller scales (in the spectral space, see Fig. 13).
Most notably, the main restriction of the eSQG framework (using SSH as input) for the density reconstruction is that both SSH and SSD have to be in phase (Lapeyre and Klein 2006; Isern-Fontanet et al. 2008). Accordingly, the L19 method is also subject to such restriction, but only at scales smaller than LC. The smaller-scale phase shift between SSH and SSD can degrade ρL19 in the shallower layers, as shown by the case study on 31 March (red line in Fig. 9c). However, such degradation mainly occurs in the mixed layer (upper 100 m), and can be ignored with pattern correlations between ρL19 and ρOFES still exceeding 0.97. Influences of phase shift at smaller scales decay rapidly with depth, so ρL19 becomes almost the same as ρW13 (black line) below the mixed layer. Statistical results (Fig. 9d) further demonstrate that the L19 method (red line) can perform as successfully as the W13 method (black line) in the upper 1000 m. The ρL19 is only slightly degraded in the first 50 m, where the pattern correlations are still quite high with values exceeding 0.93.
c. Comparison with the eSQG reconstruction
Using the same 1/30° OFES simulation, Qiu et al. (2016) have already demonstrated that the eSQG method can successfully reconstruct the subsurface vorticity field from SSHA. However, they have not analyzed the density field. Here we further evaluate the L19 method against the eSQG method in terms of both the flow and density reconstructions. Notice that we use SSHA to conduct the eSQG velocity, vorticity ζeSQG, and density ρeSQG reconstructions, according to Qiu et al. [2016, see their Eqs. (1)–(3)]. For the flow fields at the 500-m depth, the L19 method is slightly less successful: the L19-reconstructed velocities (blue arrows in Fig. 4b) is slightly less correlated with the OFES simulation (red arrows) than the eSQG solution (blue arrows in Fig. 4c), with correlations for the eSQG reconstruction reaching 0.96 for both velocity components (Table 1); pattern correlation for ζL19 (0.75, Fig. 5c) is also lower than that for ζeSQG (0.82, Fig. 14a). However, for the 500-m density reconstruction, ρL19 (Fig. 11c) is better than ρeSQG (Fig. 14b), with the pattern correlation between ρeSQG and ρOFES being slightly lower (0.9; see Table 1).
For the case on 31 March 2001, although ζL19 is not as good as ζeSQG in the 450–1000-m layers, ρL19 is better than ρeSQG in the upper 1000 m (see red and blue lines in Figs. 9a,c). Statistically speaking (Figs. 9b,d), the L19 method (red line) is quite effective: it can perform as successfully as the eSQG method (blue line) in terms of ζ reconstruction, particularly in the upper 600-m layers; it also can simultaneously achieve a better ρ reconstruction throughout the 1000-m upper ocean, especially in the shallow layers. Note that for the ζ reconstruction at 0 m (Figs. 9a,b), all the three methods present the same pattern correlation (less than 1.0 due to the existence of ageostrophic flows), because they deduce the surface ζ from the same OFES SSHA field.
5. Summary and discussion
In this paper, we employ the 1/30°-resolution OFES simulation to further evaluate the W13 method for reconstructing the ocean interior from the surface information. Results elucidate the limitation of this method: it cannot properly diagnose smaller-scale circulation structures, when the horizontal resolution used is high enough. In preparation for the upcoming SWOT satellite mission (to be launched in 2021), we further propose an extended method L19 to circumvent that limitation. Evaluation of the L19 method indicates that it can be quite satisfactorily applied to the high-resolution surface data.
The W13 limitation rests upon the deficiency of only two gravest vertical normal modes for the downward projection of high-resolution SSH: smaller scales in the (residual) SSHA, which notably contain higher-order modes and should be surface-intensified, are falsely overprojected (by only BT and BC1) into much deeper layers as spurious signals. The W13 method noticeably manifests its limitation in the ζ reconstruction below the mixed layer (especially at deeper levels), because the second-order horizontal derivatives can highlight the spurious smaller horizontal scales in the reconstructed Ψ. However, the W13 limitation is inconspicuous for the ρ reconstruction, since those spurious smaller-scale Ψ decaying slowly in the vertical can be greatly eliminated by the vertical derivative.
The L19 method retains the W13 two-gravest-mode scheme to project (residual) SSHA at larger scales (>LC, LC is set as 150 km in this study), since this scheme is reasonable for the interior PV forced by the traditional Phillips-like instability that influences larger-scale dynamics (Klein et al. 2008). Meanwhile, the L19 method exploits the eSQG framework to project SSHA at smaller scales (≤LC), because at these scales this framework can reasonably capture higher-order modes, and can properly depict the surface-intensified features that do not pertain to the traditional Phillips model (Rocha et al. 2013; Sasaki et al. 2014; Qiu et al. 2016). Statistical results over two years demonstrate that the L19 method can overcome the W13 limitation: the L19 method performs as successfully as the eSQG method (using SSH as input) in terms of ζ reconstruction in the upper 1000 m.
The advantage of the L19 method over the eSQG method using SSH as input is that the former can realize a better ρ reconstruction, especially in the mixed layer. This is because that, at larger scales (>LC), the L19 method is not subject to the eSQG restriction that both SSH and SSD have to be in phase. Although the smaller-scale (≤LC) SSH–SSD phase shift would degrade the L19-reconstructed ρL19, results demonstrate that such degradation (largely confined in the mixed layer) is quite negligible and ρL19 is still very successful throughout the upper 1000-m ocean. It is important to notice that in this study, we choose to use SSH as the surface condition for the eSQG solution [see Eq. (3)]. As suggested by Lapeyre and Klein (2006), SSH and SSD are not independent in the eSQG framework. Therefore, if SSD is used instead as the surface condition, the eSQG solution can still be fitted into the isQG framework with a modified stratification (Isern-Fontanet et al. 2008). In this case, discrepancies exist (at scales smaller than LC) between the observed total SSH and the L19-SSH solution.
Besides the upstream Kuroshio Extension subdomain (dashed box in Fig. 1), we have also conducted the experiment in another subdomain (dashed box in Fig. 15a) of the North Pacific Current. Statistical results over two years also demonstrate the effectiveness of the L19 method for reconstructing both the ζ and ρ fields in the upper 1000-m layer (see red lines in Figs. 15b,c): compared to the OFES model output, the L19-reconstructed ζ and ρ fields reach a correlation of 0.65–0.9 and 0.8–0.95, respectively. Therefore, based on comparison with model fields, the proposed L19 method presents great potential for the reconstruction of ocean interior from high-resolution surface information, including possibly the upcoming SWOT SSH. It is important to note that, the area-averaged N2(z) profile used in the L19 method may not be sufficient for the local reconstruction, especially in the upstream Kuroshio Extension region where the stratification presents a noticeable horizontal variation. One may then consider employing a local N2(x, y, z) to improve the L19 performance.
There are three points worth remarking on in concluding our work.
As suggested by several studies (e.g., Wang et al. 2018; Morrow et al. 2019), the internal gravity wave signal (that will be present in SWOT SSH) probably dominates SSH at wavelengths shorter than about 100 km. However, the OFES simulation serving as the test bed of our explorations does not include tidal forcing. Therefore, we cannot examine the extent to which the internal gravity wave can impact the L19 reconstruction (of the balanced fields). Future works that can quantify the internal tide and wave effects are needed.
At present, we tentatively set a cutoff scale LC to be 150 km, and assume that the residual SSHA at scales smaller than LC contains significantly higher-order modes. In future studies, more investigations are required for a better choice of LC.
At wavelengths shorter than LC, the L19 method (i.e., the eSQG framework) still is subject to the restriction that both SSH and SSD have to be aligned. Moreover, at submesoscales (<~50 km), the L19 method (the eSQG framework) cannot reconstruct the features that are mostly related to mixed layer instabilities (Sasaki et al. 2014), although it can represent the submesoscales resulting from frontogenesis trigged by the convergences that is induced by mesoscale eddies (Lapeyre and Klein 2006; Klein et al. 2009). New methods should be pursued in future works. A prospective way may be to relax the QG approximation using the Semi-Geostrophic approximation (see Badin 2013; Ragone and Badin 2016; Lapeyre 2017, and references therein).
We thank Shiqiu Peng, Rui Xin Huang, Bo Qiu, and Yu-Kun Qian for valuable discussions. We also thank both reviewers and the editor for their insightful and constructive comments. This work was supported by the National Natural Science Foundation of China (Grants 41806036, 41775085, and 41806035) and by Innovation Academy of South China Sea Ecology and Environmental Engineering, Chinese Academy of Sciences (ISEE2018PY05). The 1/30° OFES simulation was conducted by using the Earth Simulator under support of JAMSTEC. We gratefully acknowledge the use of high performance computing clusters at the South China Sea Institute of Oceanology, Chinese Academy of Sciences.
The W13 and eSQG Methods
In the QG theory, the streamfunction and PV are related as
where f0 is the Coriolis parameter, N(z) the Brunt–Väisälä frequency, and q the QG PV anomaly (Pedlosky 1987). The term Ψ = p(ρ0f0)−1 is the geostrophic streamfunction, with p and ρ0 respectively being the pressure anomaly and the reference density. Since the QG PV equation is linear, the total solution Ψ admits a superposition of a homogeneous solution termed the SQG solution Ψsur and a particular solution called the interior solution Ψint: Ψ = Ψsur + Ψint (Charney 1971; Hoskins 1975; Ferrari and Wunsch 2010). The term Ψsur is forced by buoyancy anomalies on the boundaries and involves zero interior PV anomaly, while Ψint is forced by nonzero interior PV anomaly and involves zero surface buoyancy, such that we have
Here b = −gρ/ρ0 is the buoyancy anomaly (with g and ρ being the gravity constant and density anomaly, respectively) whose value is neglected at the bottom boundary (H = 4000 m in this study). Equation (A2) resembles a reduction of the Eady model to the surface, and Eq. (A3) is associated with the classical QG model (such as the Phillips model). The combination of Eqs. (A2) and (A3) can be thought of in terms of the Charney-like instability (Lapeyre and Klein 2006).
a. The W13 method (Wang et al. 2013)
This method determines Ψint through a different scheme based on the vertical normal modes. Specifically, the W13 method expands Ψint in terms of these modes that are orthonormal and form a complete basis, and then truncates Ψint to the two gravest (BT and BC1) modes:
where the hat denotes the horizontal Fourier transform, (k, l) the horizontal wavenumbers. The term Fn is the traditional flat bottom normal mode (Pedlosky 1987) with An being the modal coefficient. The key point of W13 is that both Ψsur and Ψint contribute to the sea surface pressure (SSH) and bottom pressure such that
b. The eSQG method (Lapeyre and Klein 2006)
Considering an infinitely deep ocean where buoyancy anomalies vanish at great depths, The eSQG method employs a constant stratification N0 and replaces the bottom boundary conditions in Eqs. (A2) and (A3) with
to analytically solve the SQG solution Ψsur, as
where κ = (k2 + l2)1/2 is the modulus of the wavenumber vector.
Based on the argument that the interior PV and SSD anomalies are well correlated (LaCasce and Mahadevan 2006), the eSQG method approximates the total solution Ψ using an effective SQG solution obtained from the replacement of N0 in the denominator of Eq. (A8) by an empirically determined value: