## Abstract

Proposals from multiple nations to deploy air–sea flux moorings in the Southern Ocean have raised the question of how to optimize the placement of these moorings in order to maximize their utility, both as contributors to the network of observations assimilated in numerical weather prediction and also as a means to study a broad range of processes driving air–sea fluxes. This study, developed as a contribution to the Southern Ocean Observing System (SOOS), proposes criteria that can be used to determine mooring siting to obtain best estimates of net air–sea heat flux (*Q*_{net}). Flux moorings are envisioned as one component of a multiplatform observing system, providing valuable in situ point time series measurements to be used alongside satellite data and observations from autonomous platforms and ships. Assimilating models (e.g., numerical weather prediction and reanalysis products) then offer the ability to synthesize the observing system and map properties between observations. This paper develops a framework for designing mooring array configurations to maximize the independence and utility of observations. As a test case, within the meridional band from 35° to 65°S we select eight mooring sites optimized to explain the largest fraction of the total variance (and thus to ensure the least variance of residual components) in the area south of 20°S. Results yield different optimal mooring sites for low-frequency interannual heat fluxes compared with higher-frequency subseasonal fluxes. With eight moorings, we could explain a maximum of 24.6% of high-frequency *Q*_{net} variability or 44.7% of low-frequency *Q*_{net} variability.

## 1. Introduction

The Southern Ocean serves as a gateway between the atmosphere and the middepth ocean, both because its steeply sloped isopycnals bring intermediate water to the ocean surface (e.g., Marshall and Speer 2012) and because winter mode water formation mixes recently ventilated water into the ocean interior (e.g., Hanawa and Talley 2001; Cerovečki et al. 2013). The region is responsible for much of the global ocean uptake of CO_{2} (Caldeira and Duffy 2000; Sabine et al. 2004) and heat (e.g., Roemmich et al. 2015) from the atmosphere. Climate model evidence attributes much of this uptake to air–sea exchanges within the Southern Ocean (e.g., Swart et al. 2018). Large-scale net air–sea heat fluxes are effectively evaluated through measurements of ocean heat content (Roemmich et al. 2012, 2015), but direct estimates of air–sea fluxes are required to probe the mechanisms governing ocean heat uptake (Bourassa et al. 2013; Swart et al. 2019). However, due to the paucity of in situ flux observations, which are challenging to collect because of high winds and high sea state, air–sea fluxes are not well observed in the Southern Ocean, and associated reanalysis data have considerable uncertainties (Bourassa et al. 2013; Gille et al. 2016; Potter et al. 2018; Swart et al. 2019).

The Southern Ocean Observing System (SOOS) was established in 2011, with the mission of facilitating the design and implementation of a comprehensive and multidisciplinary observing system for the Southern Ocean through international collaboration (e.g., Rintoul et al. 2012; Newman et al. 2019). SOOS recognizes that observing system simulation experiments (OSSEs) are important tools to test the effectiveness of a new observing system (e.g., Errico et al. 2013). Southern Ocean OSSEs have been carried out, but they have been limited to determining the optimal number of autonomous profiling floats or air–sea exchange floats (e.g., Kamenkovich et al. 2017; Mazloff et al. 2018). Therefore, SOOS continues to advocate for the development of tools for designing an optimal observing system to measure essential ocean variables.

The Southern Ocean Observing System has identified fluxes as a priority observation gap. This led in 2015 to the establishment of the SOOS Working Group on Air–Sea Fluxes (SOFLUX) (Gille et al. 2016; Newman et al. 2019). SOFLUX has identified priorities for the coming decade, as articulated by Swart et al. (2019). These include collection of high-quality in situ point measurements necessary to understand small-scale flux variability and observation of fluxes in ice-covered regions. In this context, a number of nations, including China, India, and Brazil, have begun discussions aimed at deploying air–sea flux moorings in the Southern Ocean that would extend the geographic range of two Southern Ocean flux moorings: the Australian Southern Ocean Flux Station mooring (Fig. 1; Schulz et al. 2012), and the U.S. Ocean Observatories Initiative Southern Ocean mooring (Fig. 1; Ogle et al. 2018), which had been extended with U.K. support but was removed in January 2020. The most extensive of the new proposals is China’s concept plan for the “Big Ring,” which would space a network of six to eight moorings throughout the Southern Ocean (Chen 2018). These offer the prospect of ground truthing satellite estimates of air–sea fluxes and providing new insights into air–sea flux processes. A leading-order question is determining where these moorings should be deployed in order to best measure spatiotemporal variability in air–sea fluxes.

There are multiple criteria that one could use to design a mooring array. For example, Mazloff et al. (2018) suggested obtaining one measurement per decorrelation length scale, and Cronin et al. (2019) pursued a similar strategy. Mazloff et al. (2018) concluded that approximately 100 optimally spaced measurement platforms (i.e., one platform every 20° longitude × 6° latitude) were required to constrain the carbon and heat inventory between 35° and 70°S on time scales longer than 90 days. However, financial costs and logistical requirements would prohibit deploying and maintaining 100 uniformly spaced moorings in the Southern Ocean.

An alternative approach is to use a small number of moorings to constrain fluxes correlated with key modes of variability in the Southern Ocean, such as the southern annular mode (Marshall 2003), the Antarctic dipole mode associated with El Niño–Southern Oscillation on interannual time scales (Chen and Yuan 2004; Yuan 2004), and the wavenumber-3–4 patterns associated with synoptic-scale blocking in the Southern Hemisphere (Trenberth and Mo 1985; Cai et al. 1999; Raphael 2004; Liu et al. 2011; Manhique et al. 2011). Because these modes define large coherent patterns, they can in principle be constrained by a small number of moorings strategically positioned based on the spatial patterns of the modes. This smaller number of required observation sites is feasible with international coordination. If the true air–sea fluxes actually behaved like those modes, the spatial pattern of each mode would be fixed in time, and we would only expect to need one mooring per empirical orthogonal function (EOF) mode to capture the temporal variability of the mode. In reality, in a time-evolving system, although flux variability is dominated by high-variance modes, it is also influenced by small-variance modes. While we expect to constrain the leading modes with a small number of moorings, observations at these mooring sites might be contaminated by signals from small-variance modes, which might result in an imperfect representation of the leading modes. While we would expect to obtain a better representation of the variability by placing each mooring at an extremum of a mode, each spatial mode could have multiple extrema, and care must be taken in identifying the optimal mooring sites. A good deployment strategy will need to prioritize a sequence of mooring deployments, as we might not be able to place all moorings at one time. Also, the strategy should be adaptable to take into account the other nonmoored components of the observing system, such as autonomous surface vehicles and buoys. The objectives of this study are therefore to explore quantifiable strategies to optimize the placement of future Southern Ocean moorings, with goals of constraining a high fraction of total variance and/or measuring in regions of high local variance.

This work is a contribution to SOOS and is specifically aimed at facilitating coordination of a multiplatform observing system that will potentially also include remote sensing systems and autonomous surface platforms, such as Saildrones or Wave Gliders (Cronin et al. 2019; Swart et al. 2019). In this study, we focus on configuring a heat flux observation array, as the air–sea heat flux in the Southern Ocean is essential to the Subantarctic Mode Water formation and its interannual variability, and the oceanic heat content variability (Swart et al. 2018; Tamsitt et al. 2020). This array design approach can also be applied to constrain the carbon flux and the buoyancy flux, the latter of which depends on both the heat flux and freshwater flux (Cerovečki et al. 2011; see also section 2). We seek to design an observing system capable of detecting air–sea heat flux deviations from the mean, annual, and semiannual cycles (i.e., 12- and 6-month harmonics). As the proposed number of moorings is small and the mooring siting strategy uses correlation over the entire Southern Ocean, for this analysis we focus on large-scale (>100 km) variability and neglect episodic heat-flux events (Schulz et al. 2012; Ogle et al. 2018; Tamsitt et al. 2020), which are often localized and would require more spatially dense observations.

Details of data and methodology are given in section 2. Results are presented in three sections: first, the characteristics of the leading air–sea heat flux variability modes and their implications for mooring selection (section 3); second, the prioritization of mooring sites (section 4); third, how well the air–sea heat flux variability in the Southern Ocean can be constrained from a small number of moorings (section 5). Major findings are summarized in section 6.

## 2. Data and methods

### a. Data

The air–sea flux fields used in this study include the Japanese 55-year Reanalysis (JRA-55) data provided by Japanese Meteorological Agency (JMA) (Kobayashi et al. 2015; Harada et al. 2016) and the European Centre for Medium-Range Weather Forecasts interim reanalysis (ERA-Interim) data (Dee et al. 2011). Both JRA and ERA are widely recognized comprehensive reanalysis that cover the last half century and are actively updated. JRA is reported on a 640 × 320 grid that has a 0.56° longitude × latitude resolution, and ERA is reported on a 512 × 256 grid that has a 0.7° longitude × latitude resolution. This study uses 1979–2016 monthly averaged air–sea heat flux (*Q*_{net}, defined with positive values indicating downward heat flux into the ocean) and freshwater flux [evaporation minus precipitation (EmP)] data, which are representative of the large-scale air–sea flux variability. The JRA-55 data are available from https://rda.ucar.edu/datasets/ds628.1/, and the ERA-Interim data are available from https://www.ecmwf.int/en/forecasts/datasets/archive-datasets/reanalysis-datasets/era-interim. In this paper, we employ JRA for the primary analysis and use ERA to test the sensitivity of mooring sites. EmP data are used to calculate the buoyancy flux to assess how effectively *Q*_{net} variability can represent buoyancy flux variability.

We also take into consideration the locations of the U.S. Ocean Observatories Initiative (OOI) Southern Ocean mooring (Ogle et al. 2018) and the Australian Integrated Marine Observing System Southern Ocean Flux Station (IMOS SOFS) mooring (Schulz et al. 2012) in order to judge future mooring deployment, and we use *Q*_{net} observations at these two sites to validate the reanalyses. The OOI mooring site is located southwest of the South American continent at 54.47°S, 89.28°W, and the SOFS mooring site is located southwest of Tasmania at 47°S, 142°E. At these two sites *Q*_{net} is calculated from the 1-min meteorological observations using the COARE 3.5 flux algorithm (Edson et al. 2013), and then averaged to yield monthly terms. Mooring deployments and data processing are detailed by Schulz et al. (2012), Ogle et al. (2018), and Tamsitt et al. (2020).

A simple comparison of the monthly mean JRA and ERA data with the mooring observations (Fig. 2) suggests consistency between the reanalysis data and the mooring observations during overlapped periods. Some comparison of daily mean fluxes from reanalysis and moorings yielded similar results (Ogle et al. 2018; Tamsitt et al. 2020). As this study does not aim at evaluating the reanalysis, we do not compare heat fluxes from moorings with those from other reanalysis or on a wider range of time scales. JRA and ERA overlap with the OOI mooring observations in March 2015, January–September 2016, and December 2016, while they overlap with the SOFS mooring observations from April 2010 to February 2011, December 2011 to December 2012, May to September 2013, and April 2015 to March 2016. The difference between ERA and SOFS net heat flux observations (ERA minus SOFS) is −10.57 ± 15.66 W m^{−2} [mean ± standard deviation (SD)], and the difference between JRA and SOFS (JRA minus SOFS) is −27.39 ± 20.48 W m^{−2}. Comparison of reanalyses with OOI observations shows that the difference between ERA and OOI observation (ERA minus OOI) is −14.59 ± 6.49 W m^{−2}, and the difference between JRA and OOI (JRA minus OOI) is −25.70 ± 7.91 W m^{−2}. This indicates a negative bias in reanalysis data with respect to in situ mooring observations, with the bias in ERA smaller than in JRA. The differences of both JRA and ERA fields relative to the OOI observations have smaller standard deviations than differences of both fields relative to the SOFS observations. This could be explained by the fact that there are stronger episodic turbulent heat loss events at the SOFS site than at the OOI site (Schulz et al. 2012; Ogle et al. 2018; Tamsitt et al. 2020). As the reanalyses are only compared to observations at two mooring sites, we do not evaluate whether there are globally negative biases in the reanalysis or whether ERA is more accurate than JRA. In the following analysis, the JRA data are employed, because they have higher spatial resolution than ERA.

Exclusive economic zones (EEZ) boundaries are used to avoid siting moorings in national waters. Generally, a state’s EEZ extends 200 nautical miles (~370 km) out from its coast, except where resulting points would be closer to another country. The EEZ data (version 10) are downloaded from http://www.marineregions.org/. This dataset combines the boundaries of the world countries and the exclusive economic zones of the world (Flanders Marine Institute 2018).

The southern annular mode (SAM) and the Amundsen Sea Low (ASL) indices provide a measure of large-scale climate patterns in the Southern Ocean. The SAM, also known as the Antarctic Oscillation (AAO), describes the north–south movement and strength of the westerly wind belt that circles Antarctica, dominating the middle to higher latitudes of the Southern Hemisphere (Marshall 2003). The monthly mean AAO index since January 1979 is downloaded from http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/aao/aao.shtml. The ASL is a climatological low pressure center located over the extreme southern Pacific Ocean, off the coast of West Antarctica. Atmospheric variability in this region is larger than anywhere else in the Southern Hemisphere and exhibits significant correlations with both the SAM and El Niño–Southern Oscillation. A set of ASL indices is derived from ERA-Interim (Hosking et al. 2016), and monthly values of ASL actual central pressure are downloaded from https://legacy.bas.ac.uk/data/absl/.

### b. Data preprocessing

We preprocess our time series of monthly mean heat fluxes by first using a least squares fit to remove annual and semiannual harmonic cycles (Fig. 3a). By removing the annual cycle, we effectively assume that the annual cycle of air–sea fluxes can be well estimated based on prior knowledge of the annual cycle of incoming solar energy (e.g., Ogle et al. 2018), and that the priority in designing a flux array should be to measure the other components of the signal. Once the annual and semiannual cycles have been removed, the residual is dominated by high-frequency variability (*f* > 1.0 cpy) and displays abrupt monthly variations. We refer to this as the high-frequency signal (although it also retains low-frequency variability). The high-frequency signal is the baseline case for this analysis, and is the case considered if no frequency range is specified.

We separate the low-frequency *Q*_{net} from the full residual by applying a 19-point (19-month) Hanning filter (Fig. 3b), which suppresses variability on time periods shorter than the annual cycle (i.e., *f* > 1.0 cpy) (Fig. 3c). In the following analyses, the derived low-frequency *Q*_{net} is also spatially smoothed using an 11/cos(latitude) × 11 point box-average filter, in order to avoid the influence of some localized point maxima on mooring siting. Throughout the text, unless otherwise noted, we use the terms “low frequency” and “high frequency” to distinguish between filtered and unfiltered *Q*_{net}.

The SDs of high-frequency and low-frequency components of *Q*_{net} are mapped in Figs. 4a and 4b. Generally, the largest SD of the high-frequency component is about 50 W m^{−2}, and the largest SD of the low-frequency component is about 10 W m^{−2}.

For some applications, the critical dynamical questions focus on processes that change water density, and therefore buoyancy flux (*B*_{net}) is more important than heat flux; *B*_{net} is the sum of air–sea heat flux and freshwater (EmP) heat-equivalent flux (Cerovečki et al. 2011; Snow et al. 2016), calculated as

where *ρ*_{0} is a reference density, *c*_{p} is the specific heat for seawater, *S*_{0} is ocean surface salinity, *α* is the thermal expansion coefficient, and *β* is saline contraction coefficient. Here *α* and *β* are calculated using monthly climatological temperature and salinity data from the Comprehensive Ocean Atmosphere Dataset (Diaz et al. 2002). SDs of both high-frequency and low-frequency *B*_{net} (Figs. 4c,d) are similar to those of *Q*_{net} (Figs. 4a,b) with negligible difference (Figs. 4e,f), suggesting that in this latitude range, conclusions derived from the *Q*_{net} analysis can likely be extended to *B*_{net}, and we therefore carry out detailed analysis on *Q*_{net} only. Since salinity dominates the buoyancy budget in cold water, making temperature less important, the agreement between *B*_{net} and *Q*_{net} variations at high latitudes is somewhat surprising. We identified three possible reasons for the agreement. First, the buoyancy flux calculation does not account for the effects of sea ice redistribution, which have been shown to be important to the total buoyancy fluxes in the seasonal sea ice zone (e.g., Abernathey et al. 2016; Pellichero et al. 2018; Cerovečki et al. 2019). Second, annual and semiannual harmonic components of EmP that can contribute to *B*_{net} variability, are not considered here. Last, the magnitude of high-latitude EmP variability in reanalyses, which could influence the magnitude of *B*_{net} variability [Eq. (1)], has not yet been validated.

### c. Methods

#### 1) Empirical orthogonal functions

EOFs (e.g., North 1984; Kelly 1988; Preisendorfer 1988) are employed to identify the leading modes of variability of *Q*_{net} in the Southern Ocean. EOFs decompose a space–time field *V*(**r**, *t*), with zero mean in time, into $V\u2061(r,t)=\u2211\u2009\alpha i\u2061(t)\Phi i\u2061(r)$, where *α*_{i} represents the temporal components and Φ_{i} the spatial components of the *i*th EOF mode. Each spatial mode Φ_{i}(**r**) is orthogonal to the other modes, and each temporal component *α*_{i}(*t*) is uncorrelated to the other temporal components. The EOF modes are sorted according to their relative importance in explaining the total variance of *V*(**r**, *t*).

Because the flux data are reported on a regular longitude-by-latitude grid, there are more data values per unit area at high latitudes, which could artificially increase the variance at high latitudes. To avoid this problem, the data are weighted by multiplying observed anomalies by a weight, (cos*θ*)^{1/2} (where *θ* denotes the latitude), prior to carrying out the EOF decomposition, and the derived EOF spatial modes are then divided by this term.

#### 2) Using maximum explained variance to site moorings

EOF spatial modes could be used directly to site moorings. For example, to constrain a specific EOF mode, moorings could be placed at the extrema of the EOF’s spatial mode. However, this would result in multiple moorings expected to produce highly correlated time series.

Here, we explore an alternative approach for siting moorings. This strategy seeks locations where the flux records are formally as uncorrelated as possible in order to maximize the fraction of variance that can be explained with as few moorings as possible. We can think of the data *V*(**r**, *t*) as a matrix with dimensions *M* in space and *L* in time, with zero mean in the time domain. The projection of the data at one location **r**_{0} onto the full dataset can be computed. We denote $||V\u2061(r0,t)||2=\u2211j=1LV\u2061(r0,tj)2$ and use angle brackets for the inner product: $\u2329V\u2061(r1,t),V\u2061(r2,t)\u232a=\u2211j=1LV\u2061(r1,tj)V\u2061(r2,tj)$. The time series *V*(**r**_{0,}*t*) can be normalized to have a unit norm via the equation **e**(**r**_{0,}*t*) = *V*(**r**_{0,}*t*)/||*V*(**r**_{0,}*t*)||_{2}. Then *V*(**r**_{0,}*t*) = ⟨*V*(**r**_{0,}*t*), **e**(**r**_{0,}*t*)⟩**e**(**r**_{0,}*t*) and its variance is ⟨*V*(**r**_{0,}*t*), **e**(**r**_{0,}*t*)⟩^{2}/*L*. For the time series *V*(**r**_{i,}*t*) at a position **r**_{i}, the component correlated to *V*(**r**_{0,}*t*) is ⟨*V*(**r**_{i,}*t*), **e**(**r**_{0,}*t*)⟩**e**(**r**_{0,}*t*), and the variance of this component is Ψ_{i} = ⟨*V*(**r**_{i,}*t*), **e**(**r**_{0,}*t*)⟩^{2}/*L*. For all *i*, with *θ*_{i} denoting the latitude at position **r**_{i}, and *ω*_{i} = cos(*θ*_{i}), the area-average variance of the time series at **r**_{i} correlated with the unit vector **e**(**r**_{0,}*t*) is therefore

We seek a position **r**_{0} such that the normalized time series **e**(**r**_{0}, *t*) has the largest overall covariance with *V*(**r**, *t*) and therefore explains the largest possible fraction of the variance in the total field [e.g., minimizes the mapping error variance as shown by Eq. (10) of Bretherton et al. 1976].

Once the position **r**_{0} is identified, the projection of the time series at **r**_{0} into the full data record is used to determine the regression coefficient *γ*:

Then, the *V*(**r**_{0}, *t*)-correlated component is removed at each location **r**_{i}: *V*′(**r**_{i}, *t*) = *V*(**r**_{i}, *t*) − *γ*_{i}*V*(**r**_{0}, *t*).

This procedure can be repeated, by applying the methodology to *V*′, seeking a geographic location **r**_{1} where *V*′(**r**_{1}, *t*) is able to explain the maximum overall variance in *V*′. By repeating this process iteratively, we can select an arbitrary number of potential mooring locations with roughly independent time series. This stepwise selection of an array of *N* optimal mooring sites determines the sequence of mooring deployment and is computationally effective, although the final selected array is not necessarily optimal (i.e., among all possible choices of *N* mooring sites, the linear regression model iteratively built up from these selected sites to fit *Q*_{net} in the Southern Ocean does not necessarily have the least root-mean-square error). Less optimal solutions might occur in cases where measurement noise (or unresolved components of the signal) hinder the separation of modes of variability represented by different moorings.

#### 3) Maximum covariance analysis

Maximum covariance analysis (MCA) (von Storch and Zwiers 1999) is employed to assess how well the air–sea heat-flux variability modes in the Southern Ocean can be constrained by the *N* selected moorings. The MCA uses singular value decomposition (SVD) (Golub and van Loan 1989), and therefore in climate research it is sometimes referred to as the SVD method (Bretherton et al. 1992; Lopez et al. 2017).

We define the covariance matrix $C$_{M×N} for an *M*-point gridded *Q*_{net} field *V*(**r**, *t*) in the Southern Ocean with respect to the *Q*_{net} field *S*(**r**, *t*) at *N* mooring sites. Then, the covariance matrix $C$ can be decomposed into $C=\u2211\u2009k=1N\sigma kpkqkT$. Here, **p**_{k} is an orthonormal set of *N* vectors of length *M* called the left singular vectors, **q**_{k} is an orthonormal set of *N* vectors of length *N* called the right singular vectors, and *σ*_{k} is a nonnegative number called the singular value. The leading pattern **p**_{1} and **q**_{1} guarantees that the time expansion coefficients $a1\u2061(t)=p1TV\u2061(r,t)$ and $b1\u2061(t)=q1TS\u2061(r,t)$ have the largest covariance among many feasible choices of (**p**_{1} and **q**_{1}). Similarly, **p**_{k} and **q**_{k} guarantees that $ak\u2061(t)=pkTV\u2061(r,t)$ and $bk\u2061(t)=qkTS\u2061(r,t)$ have the largest covariance among many feasible choices of (**p**_{k}, **q**_{k}), where **p**_{k} is a unit vector orthogonal to **p**_{k−1}, **p**_{k−2}, …, **p**_{1}, and **q**_{k} is a unit vector orthogonal to **q**_{k−1}, …, **q**_{1}. The correlation between time expansion coefficients *a*_{k}(*t*) and *b*_{k}(*t*) assesses how well the corresponding *Q*_{net} variability mode can be predicted from observations at a small number of mooring sites. The results of MCA can be used to reconstruct one field from observations of another field (e.g., Bretherton et al. 1992; Lopez et al. 2017). At a time *t*_{0}, *V*(**r**, *t*_{0}) can be approximated from *S*(**r**, *t*_{0}) by utilizing the SVD modes via the formula $V\u2061(r,t0)\u2248\u2211k=1Nbk\u2061(t0)pk\u2329ak,bk\u232a/\u2329bk,bk\u232a$. In the extreme condition when the mooring sites collocate at the grid points, i.e., *M* = *N*, the MCA is equal to the EOF, and the spatial and temporal variability modes of *Q*_{net} are completely constrained by mooring observations. Similar to the EOF analysis, here the *Q*_{net} field is also area weighted by (cos*θ*)^{1/2} prior to SVD analysis, and the derived left and right spatial patterns are then divided by this term.

## 3. Spatiotemporal variability modes of *Q*_{net}

As a first step in this analysis, we use EOFs to assess the spatial structure of air–sea heat-flux variability and the number of independent modes needed to describe this variability efficiently. Dominant EOF modes are often prioritized in observations, and depending on the spatial structure, each EOF mode might require more than one mooring to discern it from other modes adequately. Ideally, these moorings would be deployed in the extrema of spatial modes.

The three leading EOF modes of high-frequency *Q*_{net} are shown in Fig. 5. The first mode displays a distinct wavenumber-3 pattern in the Southern Ocean, with the largest negative and positive amplitudes along the core of the Antarctic Circumpolar Current (ACC). This mode, which here explains 7.5% of the total variance, has also been identified in previous EOF analyses of Southern Ocean surface heat fluxes (Cai et al. 1999; Liu et al. 2011). In Fig. 5a, the largest negative amplitudes are located over the eastern Pacific sector, south of Africa, and south of Australia, while the largest positive amplitudes are in the western Atlantic Ocean sector, in the central Indian Ocean sector, and southeast of New Zealand, consistent with previous studies. This wavenumber-3 mode is quasi stationary (Cai et al. 1999) and is associated with the meridional component of the large-scale atmospheric circulation in southern high latitudes (Raphael 2004).

The second mode reflects a distinct wavenumber-4 pattern in the Southern Ocean (e.g., Wang and Dommenget 2016; Lin 2019), with high amplitude in the Atlantic sector. As Fig. 5b illustrates, the centers of the wavenumber-4 modes tend to be farther from Antarctica than the EOF 1 maxima are. This could occur because the high-wavenumber oscillation modes require longer zonal distances. The third mode reflects a wavenumber-3 pattern, with high amplitude in the Pacific and Indian Ocean sectors (Fig. 5c). Two EOF modes are required to capture the wavenumber-3 mode because it is a propagating pattern.

The EOF modes of low-frequency *Q*_{net} are presented in Fig. 6. The first EOF mode represents a trend, which is not linear, while the second and third modes represent interannual oscillations. EOF1 suggests that *Q*_{net} undergoes an overall decreasing trend, except in the regions south of Chile and southwest of Australia. This mode implies a decreasing heat flux (to the ocean) in regions with a positive spatial mode (red in Fig. 6a), implying a deceleration of warming in the Southern Ocean (Gille 2002; Armour et al. 2016). The second mode is the Antarctic dipole mode, which is linked to El Niño–Southern Oscillation (ENSO) (Chen and Yuan 2004; Yuan 2004). The third mode explains the variance near the Agulhas retroflection and in the central Pacific sector, and the amplitude of the spatial mode in the Pacific is offset by 30° longitude relative to the second mode. The temporal component of the third mode has no clear relationship to ENSO, the ASL, or the SAM.

The first three EOF modes explain 18.9% of the total variance of high-frequency *Q*_{net}, while the first eight capture 36.6% of the total variance. For the low-frequency component of *Q*_{net}, the first three modes explain 37.6% and the first eight modes explain 61.2% of the total variance. These fractions of variance explained provide a measure of the maximum variance that can be explained by a fixed number of independent moorings.

The spatial patterns of the EOF modes provide guidance on where moorings should be deployed in order to represent these modes of variability. If one EOF mode accounted for a large fraction of total variance, a single mooring deployed at an extremum of the EOF spatial pattern might be expected to determine the temporal variability of this mode. In reality, the variability at a single point generally consists of signals from different modes, and observations from multiple moorings are likely to be needed to reconstruct the overall signal. A single EOF mode could be reconstructed with less sensitivity to measurement error by deploying moorings at all of the extrema of the spatial pattern. For instance, the wavenumber-3 pattern of EOF1 could be constrained using six moorings at the six extrema in Fig. 5a. We used this information to identify a set of ideal moorings sites for representing the wavenumber-3 mode of high-frequency *Q*_{net}, and a second set to represent the second and the third EOF modes of low-frequency *Q*_{net}. We also indicate sites that best reflect the ASL and AAO based on the criterion of extrema of regression and/or correlation coefficients (Fig. 7). All of these mooring sites are shown in Fig. 8. This assessment concentrates candidate mooring sites near the Amundsen Sea, south of New Zealand, and in the middle of the Indian Ocean sector. In some regions, the mooring sites in Fig. 8 are in close geographic proximity, suggesting that heat fluxes in these regions represent the summation of multiple independent modes.

One evident problem in deploying moorings based on spatial modes is that it is difficult to decide the sequence of moorings, and moorings deployed at all the extrema of a specific EOF spatial mode are expected to capture highly correlated air–sea heat flux variability. In the next section, we address this by selecting candidate mooring sites with uncorrelated time series.

## 4. Mooring placement experiments

Here we select mooring sites following a criterion of maximizing the fraction of variance of high-frequency and low-frequency *Q*_{net} explained in the entire spatial domain—that is, the variability at an optimal mooring site, after being normalized to have unit norm, should have the maximum covariance with all other sites [see Eq. (2)]. Some restrictions are imposed. First, the moorings need to be deployed within the meridional band from 35° to 65°S, although the analysis domain is larger (i.e., south of 20°S). Second, the moorings need to be outside of any country’s EEZ, since deployments within EEZs would be in the purview of individual countries. We also avoid deploying moorings in shallow water regions (subjectively defined as <2000 m) as the moorings are also used to observe the deep ocean. Third, we require that the moorings be separated by a minimum distance equivalent to 15° longitude in order to avoid having a limited number of moorings concentrated in one ocean basin. This approach, which aims to optimize mooring location, is different from that of Mazloff et al. (2018), who sought to identify the spatial scale over which heat fluxes decorrelate.

As the high-frequency heat flux variability is less well constrained by the large-scale energy budget of the Earth system, it is therefore more challenging to determine than low-frequency *Q*_{net}. In addition, because the mooring observing system is expensive to build and often relies on short-term funding, it is unlikely to be maintained for decades. Therefore, the mooring observing system for high-frequency *Q*_{net} is taken as the baseline case. As such, the process of determining mooring sites to constrain high-frequency *Q*_{net} is illustrated in detail, and we focus on analyzing the covariability of high-frequency *Q*_{net} at the selected mooring sites with *Q*_{net} across the Southern Ocean.

### a. Placing moorings to constrain high-frequency Q_{net}

We now illustrate the process of siting moorings to constrain high-frequency *Q*_{net}. The first mooring M_{1} is placed at the position with the largest correlation with time series variability over the entire geographic domain (see section 2 for details of the methodology). The time series at M_{1} is projected onto all other sites. The portion of the time series at each point in the spatial domain that is correlated with M_{1} is subtracted out. This orthogonalizes the *Q*_{net} time series relative to M_{1}. Then, the second mooring M_{2} is sited following the same principle in order to explain the largest possible fraction of the remaining variance. This process is repeated iteratively in order to site subsequent moorings.

As seen in Fig. 9, the first mooring M_{1} is sited in the central Pacific sector along 120°W. This site corresponds to an extremum of the wavenumber 3 mode and also the place where *Q*_{net} correlates most strongly with the ASL and SAM indices. It is also a location where the second and third EOFs of low-frequency *Q*_{net} have high-amplitude signals (Fig. 8). The correlation map of *Q*_{net} with respect to this mooring site (right column in Fig. 9) suggests that this mooring captures a significant fraction of wavenumber 3 variability.

The second mooring M_{2} is placed close to the 90°W line, which is located between positive and negative extrema of the wavenumber-3 EOF1 mode, and at the negative extreme of the wavenumber-3 EOF3 mode (Fig. 5). This site is close to the OOI site (Ogle et al. 2018; Tamsitt et al. 2020). The correlation map indicates that this mooring also captures a wavenumber-3 mode, offset in longitude by 30° relative to the leading-order EOF mode (Fig. 5a). This is consistent with the fact that Southern Ocean anomalies propagate circumpolarly, often with wavenumber 3 (van Loon and Jenne 1972; Trenberth and Mo 1985; Cai et al. 1999; Raphael 2004; Liu et al. 2011; Manhique et al. 2011).

The first two moorings are sited in the Pacific sector. They both capture *Q*_{net} variance associated with the wavenumber-3 pattern of variability. After the projections of the M_{1} and M_{2} records onto Southern Ocean–wide *Q*_{net} are removed, variability in the Pacific and Atlantic sectors is less capable of explaining the total variance (see center column of Fig. 9). Hence the optimal position for the third mooring M_{3} is in the central Indian Ocean sector. The correlation map with M_{3} suggests that this mooring captures a wavenumber-4 mode, aligning more effectively with the spatial structure represented by the second EOF mode (Fig. 5b). The fourth mooring M_{4} is also located in the Indian Ocean sector and explains the wavenumber-4 mode. The fifth and sixth moorings are again in the Pacific sector (Fig. 9).

As evidenced in the center column of Fig. 9, increasing the number of moorings has progressively less impact on constraining heat flux variance. In addition, there are financial costs associated with maintaining additional moorings. For the purposes of this study, we have not identified criteria to judge the optimal number of mooring sites, as the criteria need to balance both the measured variance and financial cost. In this analysis, we will consider a maximum of eight moorings.

Figure 10a summarizes the positions of the first eight mooring sites identified through this successive variance maximization procedure. The final two moorings, which were not shown in Fig. 9, are included here: the seventh and eighth moorings are located in the Atlantic sector. The coordinates of these moorings are listed in Table 1.

An interesting phenomenon is that the SD of *Q*_{net} adjacent to the Antarctic continent, in the seasonal ice zone, is not reduced by projecting the spatial field onto an increasing number of candidate mooring locations. This suggests that although the SD is large near the Antarctic coast, the coastal *Q*_{net} variance is not well correlated with the open ocean, as is evident in Fig. 9. Thus, although the optimization procedure employed here works well for the open Southern Ocean, it is not able to simultaneously constrain high-variance fluxes in the seasonal sea ice zone.

The simplicity of this iterative method makes it a useful tool for exploring the sensitivity of mooring siting to the following different criteria:

We restricted the analysis domain to be south of 40°S (Fig. 10b) to examine how the change in the meridional domain influences the results. In this case, the first two mooring sites are unchanged while the third and fourth moorings shift from the Indian Ocean sector to the Pacific sector, and the positions of additional moorings are also modified. These shifts reflect the fact that the ACC is farther north in the Indian Ocean sector than in the Pacific sector, meaning that some of the Indian Ocean variance is suppressed when the domain is confined to the region south of 40°S.

We did not area weight the variance, i.e.,

*ω*_{i}= 1 in Eq. (2). This resulted in almost no changes in mooring sites, except that the seventh mooring in Fig. 10a was replaced by the eighth mooring in Fig. 10c, which is close to the Antarctic (Fig. 10c). In the case without area weighting, moorings tend to be sited close to Antarctica because the increased number of data points per unit area artificially emphasizes high-latitude variance.We subtracted out the components of

*Q*_{net}that were correlated with the OOI and SOFS moorings before identifying additional mooring sites. In this case, the first mooring site is almost unchanged with a little westward shift, the second mooring site now shifts to the original third mooring site, and the other mooring sites also change (Fig. 10d). These changes occur because the original second mooring site is close to OOI mooring site and is replaced by it.

These results also suggest that the first three mooring sites in Fig. 10a are generally robust and could be prioritized as key locations for obtaining measurements to constrain air–sea heat fluxes. The analysis shows a disproportionate number of moorings in the Pacific sector (Figs. 10a–d). That is because variability in this sector is well correlated with signals over a large geographic domain, as inferred from the middle panels of Fig. 9, and the Pacific sector itself is larger, so that mooring locations in the Pacific have greater skill in explaining variance over a large geographic region.

### b. Placing moorings to constrain low-frequency Q_{net}

Similar analyses were carried out for low-frequency *Q*_{net}, again with the objective of optimizing the positions of eight moorings. In this case, the optimal location for the first mooring M_{1} is in the Drake Passage, and the second M_{2} is located in the central Pacific sector (Fig. 10f), close to the location of the high-frequency first mooring (Fig. 10a). These two moorings capture the first and second EOF modes of low-frequency *Q*_{net} (Fig. 6). Other moorings are concentrated mainly in the Pacific and Indian Ocean sectors.

For the low-frequency *Q*_{net}, we extended our analysis by considering the same variant cases that we analyzed for the high-frequency field.

If the analysis domain is restricted to the region south of 40°S, the first mooring now is sited on the 120°W line in the central Pacific sector, as now the leading EOF mode has shifted from the trend mode to the Antarctic dipole mode. The other mooring sites are also modified (Fig. 10g). However, the first three mooring locations now agree closely with the first three for the high-frequency case (south of 40°S, Fig. 10b), suggesting that the same locations constrain both high- and low-frequency variance south of 40°S.

Without area weighting, the first four moorings are unchanged from the baseline case (Fig. 10f), while some of the other moorings are shifted to be closer to Antarctica (Fig. 10h).

If we subtract out the components of heat flux that are correlated with the time series at the OOI and SOFS mooring sites, the first mooring site that captures the leading mode is unchanged, while the other mooring sites shift (Fig. 10i). This suggests that the OOI and SOFS moorings contribute significant information on low-frequency

*Q*_{net}variability as demonstrated by Tamsitt et al. (2020).

The results shown in the first two columns of Fig. 10 are computed in order to maximize the total variance explained in the entire spatial domain. A potential shortcoming of this approach is that it could lead to moorings being sited at locations with low local variance, despite exhibiting high correlation with adjacent points. To evaluate this possibility, we tried an alternate approach that simply places moorings at the locations of largest local SD. This criterion of siting moorings at locations with the largest SD could also yield unsatisfactory results if highly localized point maxima were prioritized. Therefore, we did not apply this criterion to high-frequency *Q*_{net} but only to the low-frequency *Q*_{net}, which is also spatially smoothed using an 11/cos(latitude) × 11 point box-average filter, so that localized point maxima are avoided.

The results of deploying moorings based on the criterion of locally largest SD are presented in the third column of Fig. 10. There are three moorings in the Pacific sector (moorings 1, 5, and 6 in Fig. 10k), one mooring adjacent to southwest of Australia (mooring 7 in Fig. 10k). The other four moorings are located in western boundary current extensions, one in the Brazil–Malvinas confluence and three in the Agulhas Return Current (moorings 2, 3, 4, and 8 in Fig. 10k); these positions have large local SD mostly because of meanders in the position of the fronts. The mooring sites are stable when experiment conditions are changed, with the exception of the mooring site southwest of Australia, which is omitted if the analysis area is restricted to the region south of 40°S (Fig. 10l). The decision to area weight or not area weight the total variance does not influence the really local SD, so that the mooring sites are unchanged (Fig. 10m). Subtracting the components of heat flux that are correlated with the heat fluxes at the OOI and SOFS mooring sites has almost no influence on the positions of the seven moorings that are far from the OOI site (Fig. 10n), in contrast with cases evaluated using the total variance maximization criterion. Results from the largest SD criterion identify three mooring locations in the Pacific sector (moorings 1, 5, and 6 in Fig. 10k) that are consistent with locations derived to constrain the largest total variance of high-frequency *Q*_{net} (Fig. 10a). This means that the three mooring sites in the Pacific sector are robust, regardless of method, and might be prioritized for observation.

### c. Sensitivity test of mooring sites to air–sea flux products

We next asked how sensitive the mooring sites are to the choice of air–sea flux product by extending the methodology from the JRA *Q*_{net} to the ERA *Q*_{net} (last row in Fig. 10). For the high-frequency *Q*_{net}, the first four moorings are unchanged (see Figs. 10e,a), suggesting again that these sites are robust choices. For the low-frequency *Q*_{net}, the positions and prioritization of these moorings derived to constrain the overall variance are changed significantly (see Figs. 10j,f), suggesting that the criterion of maximizing the overall fraction of variance is not robust to the data product used in case of selecting moorings to constrain low-frequency *Q*_{net}. In contrast, the mooring positions derived for low-frequency ERA *Q*_{net} based on the locally largest SD criterion are almost unchanged from the JRA results, although the prioritization is altered (see Figs. 10o,k).

## 5. Covariability modes between *Q*_{net} in the Southern Ocean and in the selected mooring sites

We test the performance of our mooring placement by using the MCA method to calculate the covariability modes between *Q*_{net} in the Southern Ocean and *Q*_{net} at the eight mooring sites selected to constrain high-frequency variability (Fig. 10a). The spatial patterns and the corresponding time expansion coefficients of the five leading SVD modes are presented in Fig. 11. In each mode, the spatial patterns of *Q*_{net} in the entire Southern Ocean are consistent with the patterns for the eight selected sites (left column of Fig. 11), and the temporal expansion coefficients for the full Southern Ocean and the eight sites are well correlated (right column of Fig. 11). The correlation coefficients all exceed 0.85 for the five leading SVD modes, and they all exceed 0.79 for the eight SVD modes. These correlation coefficients reflect how well the corresponding air–sea heat-flux variability mode in the Southern Ocean can be predicted from observations at selected moorings. The first and fourth SVD modes represent the wavenumber-3 pattern as seen in the first and third EOFs of high-frequency *Q*_{net} (Fig. 5), and the second and third SVD modes represent the wavenumber-4 pattern. This indicates that, as expected, the leading-order EOFs in the Southern Ocean can be largely constrained by a small number of moorings. The fifth SVD mode represents a meridional dipole mode in the Pacific sector.

The results of MCA can be readily used to reconstruct *Q*_{net} in the Southern Ocean from observations at the selected mooring sites. Figure 12 shows the SDs of high-frequency *Q*_{net} from observations, reconstruction, and the difference between observations and reconstruction. The reconstruction captures the variance primarily near the mooring sites (Fig. 12b). The observed SDs at these mooring sites are not necessarily larger than SDs elsewhere; they are selected because they are well correlated across a large region of the ocean. The *Q*_{net} variability in some nearshore locations (e.g., black dots in Fig. 12c) are not constrained by these selected moorings despite their high variance. That is because they are not well correlated with open ocean variability. Given that some regions of high variance have significant ocean dynamical features, we might need to prioritize measurements of locally isolated variance.

The reconstructed heat flux field from eight moorings account for 24.6% of the total variance of area-weighted *Q*_{net} field. This variance fraction is small but reasonable, given that the EOF analysis demonstrates that a maximum of 36.6% of the variance of the area-weighted *Q*_{net} field can be explained using eight orthogonal modes. Although these mooring sites explain only a small fraction of variance, they are capable of constraining the dominant modes of variability, which, when combined with satellite observations and reanalysis products, are expected to improve the accuracy of air–sea flux estimates.

We also calculated the fraction of low-frequency *Q*_{net} variance explained by the eight mooring sites that have been chosen to maximize the fraction of high-frequency variance explained in Fig. 12a. For simplicity, the processes of covariability analysis and reconstruction are not presented. These sites are able to explain 39.3% of the total low-frequency *Q*_{net} variance. In contrast, eight mooring sites selected to explain the largest possible fraction of total variance of low-frequency *Q*_{net} (Fig. 10f) can explain 44.7% of the variance, while eight moorings selected based at sites of largest SD (Fig. 10k) can explain 41.7% of the low-frequency variance. These results suggest that the optimal mooring sites for constraining high-frequency *Q*_{net} are not the optimal sites for constraining low-frequency *Q*_{net}, and mooring sites selected to maximize the fraction of variance indeed explain more overall variance than sites selected because they correspond to high SD locations.

## 6. Concluding remarks

A methodology has been developed for selecting optimal sites for air–sea heat flux moorings in order to maximize the fraction of Southern Ocean *Q*_{net} variance measured by the mooring array. Compared to previous studies of flux observing system requirements that estimated the optimal number of measurement sites (Mazloff et al. 2018), this strategy is aimed at capturing leading-order variability with just a few moorings and at prioritizing the potential impact of candidate mooring sites. This strategy selects moorings that obtain air–sea heat flux measurements that are as uncorrelated as possible, and it addresses the problem of placing multiple moorings to observe highly correlated signals at the extrema of the EOF spatial patterns. This strategy is applied both to high-frequency subannual flux fields (i.e., after subtracting annual and semiannual harmonic components from monthly averaged data) and to low-frequency interannual fields. Because spatial correlation patterns are different for different time periods, an air–sea heat flux observing system focused on low-frequency *Q*_{net} might be configured slightly differently than one aimed at high-frequency *Q*_{net}. As demonstrated by the maximum covariance analysis, the moorings selected from this strategy can effectively constrain the dominant EOF modes of *Q*_{net} in the Southern Ocean.

Several alternative site selection criteria are tested. In most cases, the Pacific sector tends to have more mooring sites than the Atlantic or Indian Ocean sectors, because the flux variability in the Pacific sector is spatially correlated over a large region (as inferred from Fig. 9). This mooring selection criterion based on maximizing the fraction of overall variance explained by the mooring array neglects some sites near the coast and in the seasonal ice zone: although nearshore and seasonally ice-covered regions can have high variance, they are often not well correlated with open ocean variability. The mooring sites derived from different criteria provide guidance on how to design the mooring air–sea flux observing system. Existing mooring sites like OOI and SOFS need to be considered as they constrain considerable air–sea flux variability. The configuration of the mooring system depends on the domain in which total variance is constrained. The domain needs to encompass the spatial structure of major variability modes in the Southern Ocean. For an observing system that captures high-frequency *Q*_{net} variability, the first four mooring sites that constrain wavenumber-3 and -4 patterns are robust to area weighting the variance and to different reanalysis products. A few specific locations show up in many of the different configurations, and hence are fairly robust to different requirements. These locations are in the central Pacific sector (mooring 1 in Fig. 10a), southwest of Chile (mooring 2 in Fig. 10a, which is close to the OOI mooring site), south of New Zealand (mooring 6 in Fig. 10a), and in the central Indian Ocean sector (mooring 3 in Fig. 10a). In contrast, when our criterion maximizes the total variance of low-frequency *Q*_{net}, mooring sites are sensitive to the reanalysis product. This points out the importance of maintaining a long-term mooring observing system in order to reduce the uncertainties in low-frequency air–sea flux variability. A challenge raised by these results is that the optimal mooring placement differs for capturing high- and low-frequency variability. Therefore, plans to develop a network to capture a wide range of scales would need to balance these requirements.

Siting moorings at locations that maximize the overall fraction of heat-flux variance might not be optimal if there are nearly comparable large variance regions that remain unsampled and if these regions feature important ocean dynamics. To prioritize some regions of high air–sea heat flux variability, a criterion based on largest local variance is also proposed and applied to low-frequency *Q*_{net}. The mooring sites selected to constrain local variance are different from those selected to maximize the global variance, but these sites are less sensitive to reanalysis products. Moreover, siting moorings at locations of maximum local variance might be unwise, if multiple locations emerge with comparably large variances and if there are already moorings close to the identified target locations. In practice, rather than simply following a single rule of maximizing variance explained, a realistic mooring deployment strategy needs to consider the potential for constraining both the local variance and the total variance, while also considering the distribution of existing moorings.

Uncertainties in reanalyses might lead to mooring locations that are not truly optimized to observe the desired features, as is evident for the mooring observing system configured to constrain low-frequency *Q*_{net}. To update our results, newly released reanalysis data that have improved data quality should be used, such as the fifth major global reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ERA5), which has higher spatial and temporal resolution than ERA-Interim (https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5). As time evolves, the spatial pattern of trends and variability may change, and impact the significance of mooring sites optimized to observe past air–sea flux variability. Nevertheless, observations from the proposed mooring sites determined from available reanalysis data will still enrich our knowledge about the air–sea fluxes and reduce uncertainties in reanalysis. These will thus contribute to developing and optimizing future mooring observing systems.

This air–sea heat flux observing system is designed to maximize the fraction of air–sea heat-flux variability explained in the Southern Ocean. However, moored buoys are also used to observe other flux variables including carbon, momentum, and freshwater fluxes. Using the same strategies presented here, a multivariate observing system could be developed to maximize the overall variability explained for multiple variables. In this situation, each variable should be normalized to have unit area-average variance (or to appropriately weight the importance of different measured variables relative to each other), with the corresponding $\Psi \xaf$ calculated according to Eq. (2). The optimal mooring site would correspond to the position **r**_{0} that maximizes the sum of the $\Psi \xaf$ values for the different variables. The effectiveness of our proposed mooring sites could be further evaluated through an OSSE, for example, to test whether one candidate set of sites helps to reduce the error in reanalysis products more effectively than other sites. Mooring deployment planning can be further informed by accounting for other components of the observing system, including ships and autonomous platforms. In situ flux observations are important contributors to assimilation efforts for numerical weather prediction and have the potential to improve the skill in predicting subseasonal to seasonal weather. Therefore, numerical weather prediction requirements could also be used to help optimize mooring siting.

## Acknowledgments

We thank the three anonymous reviewers for their constructive comments. This study is a contribution to the Southern Ocean Observing System, by the Southern Ocean Fluxes (SOFLUX) and Observing System Design Working Groups. The analysis codes used to site moorings are available via the link https://github.com/ywei-sio/SO_Mooring_Flux_Array. The ERA-Interim data can be downloaded at https://www.ecmwf.int/en/forecasts/datasets/archive-datasets/reanalysis-datasets/era-interim. The JRA-55 data can be downloaded at https://rda.ucar.edu/datasets/ds628.1/. OOI flux data can be obtained from the NSF Ocean Observatories Initiative Data Portal (http://ooinet.oceanobservatories.org), and SOFS flux data can be obtained from the Integrated Marine Observing System (IMOS) Australian Ocean Data Network portal (https://portal.aodn.org.au). STG and MRM were supported by the National Science Foundation (NSF) Grants OCE-1658001 and PLR-1425989. DC was supported by the Chinese Academician Science and Technology Consulting Project Grant 2018-XZ-12-02 and the National Natural Science Foundation of China (NSFC) Grant 41621064. YW was supported by the NSFC Grant 41806024 and International Postdoctoral Exchange Fellowship Program (20190016) by the office of China Postdoctoral Council. LN is supported by the Australian Research Councils Antarctic Gateway Partnership (SRI40300001). VT acknowledges support from the Centre for Southern Hemisphere Oceans Research (CSHOR); CSHOR is a joint research Centre for Southern Hemisphere Ocean Research between QNLM and CSIRO. SS is supported by a Wallenberg Academy Fellowship (WAF 2015.0186) and from the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement 821001 (SO-CHIC). The authors have no conflict of interest to declare.

## REFERENCES

*Nat. Geosci.*

*Nat. Geosci.*

*Bull. Amer. Meteor. Soc.*

*J. Climate*

*Deep-Sea Res. Oceanogr. Abstr.*

*J. Climate*

*Science*

*J. Climate*

*J. Phys. Oceanogr.*

*J. Climate*

*J. Climate*

*Front. Mar. Sci.*

*Quart. J. Roy. Meteor. Soc.*

*WMO. Bull.*

*J. Phys. Oceanogr.*

*Quart. J. Roy. Meteor. Soc.*

*Science*

*Eos, Trans. Amer. Geophys. Union*

*Matrix Computations*

*Ocean Circulation and Climate: Observing and Modelling the Global Ocean*, G. Siedler, J. Church, and W. J. Gould, Eds., International Geophysics Series, Vol. 77, Academic Press

*J. Meteor. Soc. Japan*

*Geophys. Res. Lett.*

*J. Geophys. Res. Oceans*

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*J. Meteor. Soc. Japan*

*Atmosphere*

*J. Climate*

*Geophys. Res. Lett.*

*Int. J. Climatol.*

*J. Climate*

*Nat. Geosci.*

*J. Geophys. Res. Oceans*

*Front. Mar. Sci.*

*J. Atmos. Sci.*

*Geophys. Res. Lett.*

*Nat. Commun.*

*Bull. Amer. Meteor. Soc.*

*Principal Component Analysis in Meteorology and Oceanography*

*Geophys. Res. Lett.*

*Oceanography*

*Nat. Climate Change*

*Nat. Climate Change*

_{2}

*Science*

*Geophys. Res. Lett.*

*J. Climate*

*Nat. Geosci.*

*Front. Mar. Sci.*

*J. Climate*

*Mon. Wea. Rev.*

*J. Geophys. Res.*

*Statistical Analysis in Climate Research*

*Climate Dyn.*

*Antarct. Sci.*

## Footnotes

Denotes content that is immediately available upon publication as open access.

Polar2018, Davos, Switzerland, Institute for Snow and Avalanche Research