SST Dynamics at Different Scales: Evaluating the Oceanographic Model Resolution Skill to Represent SST Processes in the Southern Ocean

In this study we demonstrate the many strengths of scale analysis: we use it to evaluate the Nucleus for European Modelling of the Ocean model skill in representing sea surface temperature (SST) in the Southern Ocean by comparing three model resolutions: 1/12 ◦ , 1/4 ◦ , and 1 ◦ . We show that while 4–5 times resolution scale is sufficient for each model resolution to reproduce the magnitude of satellite Earth Observation (EO) SST spatial variability to within ± 10%, the representation of ∼ 100-km SST variability patterns is substantially (e.g., ∼ 50% at 750 km) improved by increasing model resolution from 1 ◦ to 1/12 ◦ . We also analyzed the dominant scales of the SST model input drivers (short-wave radiation, air-sea heat fluxes, wind stress components, wind stress curl, and bathymetry) variability with the purpose of determining the optimal SST model input driver resolution. The SST magnitude of variability is shown to scale with two power law regimes separated by a scaling break at ∼ 200-km scale. The analysis of the spatial and temporal scales of dominant SST driver impact helps to interpret this scaling break as a separation between two different dynamical regimes: the (relatively) fast SST dynamics below ∼ 200 km governed by eddies, fronts, Ekman upwelling, and air-sea heat exchange, while above ∼ 200 km the SST variability is dominated by long-term (seasonal and supraseasonal) modes and the SST geography.


Introduction
The choice of spatial and temporal resolution of oceanographic models is a trade-off between computational cost and the realistic representation of physical, chemical, and biological phenomena. In this study, we address several questions regarding the optimal spatial resolution of oceanographic models and their forcing data sets from the point of how the model represents the sea surface temperature (SST) in the Southern Ocean (SO) above the submesoscale (>10 km). The SST field has an obvious advantage in that there are large volumes of readily available satellite Earth Observation (EO) data, but SST is also of major interest because it provides an important ocean feedback to the atmosphere and hence weather (e.g., Byrne et al., 2015;Frenger et al., 2013). SST variability is also indicative of fronts and baroclinic processes in the upper ocean, while frontal structure is a key feature responsible for the ocean-atmosphere exchange of heat, carbon, and other trace gases. The exchange of gases between ocean and atmosphere is particularly important in the SO where water is ventilated between the deep ocean and the surface. This makes the SO a hugely important region for understanding Earth's climate and how it may change in future decades (e.g., Carter et al., 2008;Rintoul, 2018). However, the SO is chronically undersampled and oceanographic models play a key role in our understanding of the region. Therefore, in climate studies it is critically important to evaluate and improve the accuracy of how oceanographic models represent SST in the SO. This accuracy depends crucially on the model spatial resolution, as both mesoscale and submesoscale eddies play critical roles in the vertical and meridional transport of heat (e.g., Bennett & White, 1986;Bôas et al., 2015;Bishop & Bryan, 2013;Hausmann & Czaja, 2012;Qiu & Chen, 2005;Roemmich & Gilson, 2001;Su et al., 2018).
Despite a number of existing skill evaluation metrics used in oceanography (Allen & Somerfield, 2009;Doney et al., 2009;Jolliff et al., 2009;Radić & Clarke, 2011;Saux Picart et al., 2012;Skákala et al., 2016;Stow et al., 2009;Taylor, 2001), studies comparing the performance of the same model with different reso-lutions (Bricheno et al., 2014;Bryan et al., 2014;Chen et al., 2013;Hewitt et al., 2016;Winton et al., 2014) have only recently appeared. In this study we compare the results of three different resolution (1/12 • , 1/4 • , and 1 • ) runs of the Nucleus for European Modelling of the Ocean (NEMO) model against satellite EO SST. The model representation of SST largely depends on the quality and resolution of the model input forcing. We therefore analyze the optimal resolution for the model input forcing data deemed likely to have a substantial impact on SST (later called "model input SST drivers"): latent, sensible, long-wave heat fluxes (sum of these is the atmospheric heat flux, AHF), short-wave incoming radiation (SWR), zonal and meridional wind stress (ZWS, MWS), wind stress curl (WSC), and bathymetry. The impact of model input forcing can be understood through the processes that influence SST across a wide range of spatial and temporal scales. These processes and their links to the model input SST drivers are summarized in Figure 1 (see also the figure caption).
A convenient tool to capture SST processes across a range of scales is analyzing the SST variability scale dependence (SST scaling analysis, e.g., Batchelor, 1959;Obukhov, 1949;. The SST scaling analysis originates in the theory of turbulent energy spectra (e.g., Gurvich & Yaglom, 1967;Kolmogorov, 1941;Novikov & Stewart, 1964;Panchev & Leith, 1972), which can be extended to spectral representation of tracers passively advected by turbulent In this study, we focus on the range of scales above the submesoscale, by which we simply refer to our limit of spatial resolution. The EO data we use have high spatial resolution as they incorporate infrared (IR) measurements. However, owing to their dependence upon microwave SST during period of heavy cloud cover, the resolving capability of the data is degraded. We estimated the effective resolution to be approximately 15-20 km. (That is, the EO data resolve eddy-like anomalies with diameters of 15-20 km or equivalently wavelengths twice these values.) Thus, our use of the term submesoscale is not very different from the more classical definition based on nondimensional numbers-i.e., flows characterized by order 1 Rossby and Richardson numbers (McWilliams, 2016;Thomas, 2008) and which often corresponds to spatial scales ≲10 km and time scales of hours to a day (Buckingham et al., , 2017Erickson & Thompson, 2018;McWilliams, 2016;Thompson et al., 2016). This understanding of submesoscales should be distinguished from another widely used definition, which is based on the scale at which geostrophy begins to break down, approximately below the radius of deformation (Sasaki et al., 2014;Su et al., 2018;Torres et al., 2018). Due to lack of sufficient SST data at these smaller scales, we do not address scales ≲15 km. This is not to say that such scales are unimportant, as they can have a profound effect on climate. For example, a recent study  comparing two high-resolution global ocean models (at 1/24 • and 1/48 • ) concluded that increasing model resolution to ∼2 km can have a major impact on representation of vertical heat fluxes, SST, and thus climate.
The precise formulation of the three questions we aim to answer in this study is as follows: Question 1. How does model resolution impact on model skill in representing SO SST above the submesoscale? Question 2. Which model input drivers need to be well resolved in order to reproduce the observed SO SST variability at a given scale? Question 3. What does the scaling of SST variability tell us about SST dynamics above the submesoscale?
The different aspects of SST scaling analysis used in this study include the following: • Capturing how well the model reproduces both the magnitude of, and spatial patterns in, SST variability at a range of spatial scales. • Identifying how SST and model input drivers (SST drivers) spatial variability distributes over a range of spatial scales. This allows the deduction of the necessary spatial resolution of each individual model input driver in order for SST spatial variability to be correctly reproduced. • Using scaling power law approximations in order to describe spatial variability scaling. Based on this the variability scaling is split into characteristic scaling regimes, separated by "scaling breaks." Power laws, together with their scaling exponents, enable the estimation of how efficiently small scale variability in SST can be captured by models by increasing the model resolution.

Delimiting and Defining the SO
The SO spreads across the whole range of longitudes, whereas in the south it is bounded by the Antarctic continental land mass. The definition of the SO northern boundary is to some extent ambiguous, in the literature being variously defined with a fixed (often somewhere between 45 • S and 60 • S) or a variable northern boundary. The physics of the SO is largely associated with the Antarctic Circumpolar Current (ACC), which circulates around the globe in the eastward direction. The northern boundary of the ACC varies zonally, seasonally, and intra-annually (Carter et al., 2008;Moore et al., 1999;Orsi & Whitworth, 2005;Orsi et al., 1995;Whitworth, 1988). Zonally, the ACC is narrowest at the Drake passage; in the Atlantic basin it expands substantially toward the southern coast of Africa (Whitworth, 1988). To completely capture the physical impacts of the ACC, it is desirable that the region of our study envelopes the whole of the ACC and its associated fronts. Focusing only on the ACC, it would seem reasonable to study a region with a varying northern boundary following the shape of the ACC. However, physics can be strongly dependent upon latitude: for example, the Rossby radius of deformation, or incoming short-wave radiation are meridionally dependent and this can have important impacts on both SST and model input driver scaling. Using a variable northern boundary could therefore substantially increase the uncertainty in the interpretation of the scaling analyses. For this reason we decided to choose a region with a constant northern boundary by estimating the ACC width using the NEMO 1/12 • resolution model outputs for both the zonal component of wind stress ( Figure 2) and ocean surface current velocity (not shown here) for 2010. Based on these data, we placed the northern boundary for the studied region at a latitude of 30 • S. A similar choice for the SO northern boundary was made in other modeling projects, such as The SO State Estimate (e.g., Mazloff et al., 2010).
In order to focus on specific regional ocean basins, we split the SO into three sections: Atlantic, Pacific, and Indian. The regional boundaries ( Figure 2) were Pacific 152 • E to 67 • W; Atlantic 67 • W to 33 • E; and Indian 33 • E to 152 • E. All regions had latitudes between 77 • S and 30 • S. The subdivision of SO into the three main oceanic basins corresponds to important zonal differences in the ACC flow and SO physics (e.g., Tamsitt et al., 2016).

The Model Data
This paper uses the NEMO global model, run for 2010 with three near-consistent forced, ocean-only simulations at different resolutions, using a tripolar ORCA grid (ORCA12, ORCA025, and ORCA1; see Figure 3).
The 1/12 • resolution configuration is described by Marzocchi et al. (2015) with the 1/4 • and 1 • resolution configurations made as consistent as possible with the 1/12 • reference. The three forced ocean-only simulations were undertaken for the 1978-2010 period using the NEMO ocean model and LIM sea ice models (as described by Marzocchi et al., 2015). The three simulations used the same ocean forcing: DRAKKAR Forcing Set 4.1 (DFS4.1, Brodeau et al., 2010). These model configurations are reasonably similar to ocean components of the HadGEM3 coupled model described in Hewitt et al. (2016), the main difference being that they employ the LIM2 (Bouillon et al., 2009;Fichefet & Maqueda, 1997) rather than CICE (Hunke et al., 2010) sea ice model. The initial value conditions were taken from the World Ocean Atlas climatological fields (Antonov et al., 2006;Locarnini et al., 2005). All the ORCA configurations used 75 vertical layers. The thickness of the top layer is 1 m, with 22 layers in the first 100 m and a total maximum layer thickness of 250 m (Marzocchi et al., 2015). The high-resolution 1/12 • configuration employs a nonlinear free surface; free slip condition at topographic and coastal boundaries; mixed bi-Laplacian and Laplacian viscosity; and Laplacian diffusion, rotated along isopycnals (Marzocchi et al. (2015)). Nearly consistent simulations were subsequently undertaken at 1/4 • and 1 • ocean resolution, with only essential resolution-related changes summarized in Table 1. It is noteworthy that the tracer diffusion coefficients in Table 1 are typically poorly constrained (Gnanadesikan et al., 2013) but play an important role in shaping the spatial SST patterns. The vertical mixing of tracers and momentum is parameterized by a modified version of the Gaspar et al. (1990) turbulent kinetic energy (TKE) scheme (described in Madec & The NEMO team, 2015;Megann et al., 2014, with some further details in Marzocchi et al., 2015). Full description of the TKE and Gent-McWilliams (GM) eddy parameterization used in the 1 • configuration of this study can be found in Appendix A.18 of Danabasoglu et al. (2014). The 1/4 • and 1/12 • configurations used the same TKE scheme but did not employ a mesoscale parameterization. (Mesoscale is here understood as in quasi-geostrophic framework: it is a scale similar to, or larger than, the first Rossby deformation radius, where Rossby number is usually much smaller than 1; e.g., Su et al., 2014;Su & Ingersoll, 2016.) Although the runs starting in 1978 did not reach full equilibrium, the simulations were long enough for the upper ocean and mixed layer to be equilibrated (Gregory et al., 2015;Held et al., 2010). The differences in the model configuration beyond the model resolution are not expected to have a major impact on the SST field. The 1/12 • model resolves mesoscale eddies (e.g., Hallberg, 2013), 1/4 • resolution permits mesoscale eddies, and 1 • model has impact of mesoscale eddies parameterized. The higher-resolution model has better resolved frontal and boundary currents, for example the ACC (Hewitt et al., 2017). Higher-resolution also better resolves topography thereby improving model dynamics. 10.1029/2018JC014791 The model runs were performed by the National Oceanography Centre, Southampton, UK, with successive 5-day average SST outputs stored at the UK Centre for Environmental Data Analysis JASMIN servers. We extracted the SST model outputs from the upper model layer (0.5-m depth).

The EO Data
We compared the model outputs with the National Aeronautics and Space Administration (NASA) Group High Resolution Sea Surface Temperature (GHRSST) Remote Sensing Systems (REMSS) global Level 4 (L4) SST analysis (www.remss.com). For a comparable model skill analysis the EO data need to have the same resolution as the highest model resolution. The standard available satellite data are either based on (i) microwave spectra measurements, which gives good coverage, but coarser resolution (∼25 km) or (ii) IR spectra measurements, which have better resolution (∼1 km), but a large amount of missing data mostly due to the presence of clouds.
To consistently compare the model 5-day averages with the EO data, we averaged the daily EO products over a 5-day time period. Satellite products with a large number of missing data pose a nontrivial problem for 5-day averaging as the highly variable presence of clouds can generate spatially skewed 5-day distributions.
The most optimal solution is to use L4 merged IR-microwave data sets, such as the GHRSST REMSS product, which have the desirable resolution and almost full coverage. The penalty for this is that those satellite products are based on statistical techniques, such as Optimal Interpolation (OI; Reynolds & Smith, 1994;Reynolds et al., 2007), which can create artificially smooth distributions at the scale of data resolution (Dong et al., 2006). The product used here is based on an OI between the NASA Advanced Microwave Scanning Radiometer-Earth Observing System, the Tropical Rainfall Measuring Mission Microwave Imager, and IR sensors such as the MODerate Resolution Imaging Spectroradiometer. It is produced daily on a 1/11.5 • grid, and we averaged it over a 5-day time scale precisely matching the model outputs. The length-correlation radius used for the OI is roughly 100 km (http://www.remss.com/measurements/sea-surface-temperature/ -oisst-description/). The effective resolution of GHRSST REMSS is limited by the resolving capability of IR and microwave SST, with the latter becoming more important during periods of persistent and spatially extensive cloud cover. The minimal spatial scale resolved by microwave SST is on the order of 40-50 km independent of latitude (Chelton & Wentz, 2005), while the minimal spatial scale resolved by IR SST is an order of magnitude smaller-on the order of 3 km for images from the Advanced Very High Resolution Radiometer, MODerate-resolution Imaging Spectroradiometer, and Visible Infrared Imaging Radiometer Suite. (The resolving capability of IR SST is currently a subject of study [Wu et al., 2017], with uncertainty being dominated by the correction for atmospheric water vapor. Such corrections can be improved upon and, in some cases, have resulted in images with unprecedented resolution. Despite these advances, identifying the resolving capability of IR SST is nontrivial.) Since time averaging removes small scale noise, the 5-day-averaged GHRSST REMSS data are expected to have increased signal-to-noise ratio when compared to the original 1-day data set. Based on spectral methods (not shown here) we estimated the effective resolution of the 5-day GHRSST REMSS to be around 15-20 km. The quality of the EO data is lowest in winter, as increased cloud cover skews satellite data (at least south of 50 • S) toward using microwave sources with coarser spatial resolution increasing the EO SST errors. Furthermore, in winter the satellite data south of 55-60 • S are likely to be significantly contaminated by floating icebergs (the ice cover in SO is largest around August; see Holland, 2014).
The smoothing originating from OI could reduce small scale variability in the EO data. The magnitude of this effect is difficult to estimate, but as discussed later, the results seem to indicate that there is only a minor impact of this on the SST scaling analysis. Furthermore, we assume that the smoothing effect of OI is partially absorbed by the 5-day time averaging, which also smooths fields at smaller scales. Because the smoothed EO data will have less small scale dynamics, we assume that the 5-day time averaging has larger impact on the NEMO 1/12 • outputs partly leveling out the artificial differences between EO and NEMO 1/12 • SST originating from the processing of the EO data.
The 5-day averaging also removes the potentially large diurnal cycle present in the satellite data, since the satellite measures the ocean skin (<1 mm) temperature. With the diurnal cycle removed, the 5-day EO data averages can be thought as representative of the model SST 5-day averages (for some discussion on the relationship between ocean skin and ocean bulk temperature see e.g., Minnett, 2003).

Methods
The methodology used to answer questions 1-3 is fully explained in supporting information S1. Here we briefly summarize the most important concepts used throughout the paper: The magnitude of spatial variability (mosv) is for (arbitrary) field f defined as The averaging ("⟨⟩") in equation (1) runs through all 2010 outputs and includes the absolute differences in between all pairs of points separated by horizontal scale within the spatial horizontal 2-D domain labeled by coordinates ⃗ x. Following the multifractal literature (e.g., Lovejoy & Schertzer, 2013;Schertzer & Lovejoy, 1987, we call the f absolute differences "f increments." We also focus on the magnitude of time variability (motv) f defined as (2) The averaging in equation (2) takes into account all the spatial points ⃗ x and all the 2010 5-day outputs separated by the time scale (time has a 5-day resolution).
We will show that the mosv scales (i.e., depends on spatial scale) as a piecewise power law: which means that f scales with two distinct power law regimes separated by the scaling break at the scale s . The two power laws are characterized by two distinct scaling exponents H 1 and H 2 (for the details, see supporting information S1, section S1).
We suggest to place the optimal resolution of an SST model input driver on the scale N , where we resolved N% from the SST driver L scale subgrid variability (for details, see supporting information S1, section S1). For piecewise power laws (equation (3)), one can derive a simple analytical expression for the SST model input driver d field optimal resolution scale N (question 2) as a function of resolved variability threshold N:

Model Input Drivers for SST Dynamics
In this section we give a more detailed derivation of the SST input drivers and the processes controlling the SST distribution above the submesoscale. We also discuss the data used for the model input SST drivers, as well as the associated uncertainties of the analysis.
The variability of SST is controlled by the mixed layer heat budget equation that can be written following Frankignoul (1985) and Dong et al. (2007): where T is temperature in the mixed layer, t is time, h is the mixed layer depth, ⃗ v is the horizontal velocity vector and w e the entrainment velocity, Q is the total net heat flux through the ocean surface, and Q + is the heat flux through the bottom of the mixed layer. Furthermore, is water density, C p its specific heat capacity, and its horizontal diffusion coefficient. Dynamical equations such as equation (5) provide a natural connection between the SST and the SST driver mosv. For detailed arguments about how this connection relates to the SST field dynamics, see section S3 in supporting information S1.
While equation (5) represents the heat budget for the mixed layer, we are making the assumption here that SST is representative of the mixed layer temperature given the sustained upper water mixing experienced in the SO as a result of the sustained winds and energetic waves and surface currents (Gille, 2005). This assumption has been verified by the model data (see supporting information S1: Figure S5). Starting from the left, equation (5), shows that the evolution of SST is influenced by horizontal advection, diffusion, entrainment processes at the base of the mixed layer, and surface-air sea heat fluxes.
whereẑ here is the unit vertical vector and ⃗ is the wind stress vector. The Ekman meridional transport was shown by Dong et al. (2007) to play a dominant role over geostrophic advection in their SST budget analysis (see also Rintoul & England, 2002;Sallée et al., 2008); however, in the Pacific (i.e., near the western boundary) the heat transport through the geostrophic flow becomes important (Tamsitt et al., 2016). The large-scale circulation in the SO has been the topic of intense research over the last two decades (e.g., Peña-Molino et al., 2014;Tansley & Marshall, 2001;Vivier et al., 2005). A significant percentage of the circulation can be considered barotropic (Peña-Molino et al., 2014;Vivier et al., 2005), and bathymetry plays an important role in shaping the patterns of circulation as well as the transport associated with the largest SO current, the ACC (Tamsitt et al., 2016), and hence, it has been included in our analyses. Similarly, the ACC transport has been shown to be sensitive to the zonally integrated wind stress (e.g., Tansley & Marshall, 2001). Our study therefore includes the following input SST variability drivers associated with horizontal advection: ZWS, MWS and bathymetry. 2. The diffusive coefficient ( , equation (5): second term on the right side) homogenizes SST distributions via horizontal mixing, and conceptually, it can be considered proportional to the length scale E associated with turbulent eddies (e.g., Visbeck et al., 1997). The eddy driven mixing is extremely important in the SO (e.g., Frenger et al., 2015;Foppert et al., 2017;Mazloff et al., 2010;Olbers et al., 2004;Rintoul & Sokolov, 2001). Although mesoscale turbulent eddies are not a model input SST driver, one can use universal multifractal theory (see supporting information S1: section S2) and determine the scaling profile of SST magnitude of variability for a hypothetical situation in which the SST was driven purely by mesoscale turbulent eddies (as a passive tracer). This scaling profile can be used to determine turbulence induced scales of SST variability. We used the universal multifractal theory and included the SST variability associated with turbulent eddies in our study. 3. The entrainment velocity w e (equation (5): third term on the right side) is largely driven by Ekman pumping (Dong & Kelly, 2004;Dong et al., 2007) such that The Ekman pumping is of special relevance to high latitudes where the wind is strong (Frankignoul, 1985). This suggests that the WSC (which would otherwise also impact the circulation in the ACC Vivier et al., 2005) should be considered as one of the input SST drivers. WSC is therefore included in this study to represent SST variability associated with entrainment processes. 4. The net surface heat flux Q (equation (5): fourth term on the right side) is composed of the SWR, long-wave outgoing thermal radiation, latent and sensible heat flux terms. In the SO the SWR is the largest component (Dong et al., 2007;Hausmann & Czaja, 2012;Hausmann, Czaja, et al., 2016) and we have therefore included it separately as an SST input driver. The other components have been considered together (only the sum of these components was available in the data) as AHF. Q + represents the total heat flux below the mixed layer, which according to Dong et al. (2007) is small compared to the uncertainty associated with Q and therefore has been ignored in our analyses.
The impact of heat fluxes and entrainment on SST depends on the mixed layer depth h (equation (5)). The mixed layer depth in turn depends on entrainment processes, heat fluxes (thermal expansion of the water column), but especially it depends on wind mixing (third power of wind speed, or | | 3/2 ; Frankignoul, 1985). We therefore included the scaling of cubic wind speed as an SST input driver. We are not attempting to close the surface mixed layer budget such as the work done by Chi et al. (2014), Dong et al. (2007), and Moisan and Niiler (1998) but rather evaluate which of the variables that have a leading control on the budget have a measurable impact on the scale dependent resolution of SST variability. While some of our selected drivers (SWR, WSC) have a direct representation in the mixed layer heat budget and SST, others, such as bathymetry, have a wider relevance in the full model solution and their impact in SST variability is highly nonlinear and important even if it cannot be isolated within the heat budget equation.

10.1029/2018JC014791
We obtained the bathymetry data from the British Oceanographic Data Centre General Bathymetric Charts of the Oceans (GEBCO; Weatherall et al., 2015) 1/12 • 2014 data set. The 1/12 • NEMO model outputs provide data for all the selected atmospheric drivers. They originate from 1 • DFS4.1 ocean forcing set (Brodeau et al., 2010). In the simulations, the DFS data were kept on its original 1 • grid but with values over land removed (using the DFS land mask) and flood filled with oceanic values (Kara et al., 2007). These data (e.g., 6-hourly 10-m winds, daily SWR) were read in by the model and interpolated during the model run to the NEMO grid. Wind speed components were interpolated using bicubic interpolation weights and rotated where necessary to fit the local grid orientation; all other components used bilinear weights. The interpolated values provided atmospheric field values to the CORE bulk formula which, taking into account the current oceanic and sea ice state, provided surface fluxes into the ocean at each grid cell. This means the AHF obtained from model outputs depends on 1/12 • model skill to represent SST. The model outputs are fairly unique in that they provide data for all of the selected SST atmospheric drivers. Furthermore, the data are internally consistent (e.g., the long-wave radiation, latent and sensible heat fluxes were calculated from the SST forced by the winds and the short-wave radiation). However, the coarse (5-day) time resolution of the model outputs smooths lower scale variability. This is not necessarily a problem, so long as the smoothing has a similar effect on all of the SST input drivers. For example, this is not going to be the case for bathymetry, since bathymetry does not depend on time. This means the bathymetry optimal resolution scale N (equation (4)) might be larger than implied by our analysis. Even with this uncertainty, using the model outputs is still the optimal choice. Otherwise we would have to limit the number of model input SST drivers (the analysis would lose much of its relevance) or compare SST atmospheric drivers with different time resolutions. Yet doing so would only increase the uncertainty in our analysis.
The model outputs for the SST drivers have high spatial resolution (1/12 • ); however, the interpolation used by the model might have importantly influenced how the SST driver variability is represented below the 1 • scale. Interpolation does not add any new variability below the 1 • scale, and the missing variability might lead to artificially steep scaling profiles (artificially large scaling exponents H at scales below 1 • ). Some discussions on this can be found in the caption to Figure S6 in supporting information S1, where we validated wind data from model outputs by comparing wind speed to an observational 1/4 . This has relevance for ZWS and MWS as all ZWS, MWS, wind speed, and third power of wind speed (collectively called wind fields) have very similar scaling (see supporting information S1: Figure S7). As a consequence, only results for ZWS are presented in this paper.
The WSC was not a specific model output, we instead calculated it from the model outputs for wind stress components at the 1/12 • scale. To remove the potentially significant impact of noise on the wind stress gradients, we smoothed the wind stress components at the 1/12 • scale using a 40-km radius moving average prior to calculating the wind stress derivatives. We calculated that the smoothing removed about 20-30% of 8-km scale magnitude of variability from the wind stress vector components.
It is instructive to compare the optimal resolution scale N of the model forcing with the N for the eddy-induced SST variability. The effect of isotropic, but inhomogeneous turbulence (in the inertial range), on the SST mosv can be estimated from universal multifractal theory (e.g., Gagnon et al., 2006;Lovejoy & Schertzer, 2013;Schertzer & Lovejoy, 1987 applied to passive scalars. Ocean is a highly stratified fluid with horizontal dimensions much larger than the vertical dimension, so it is highly desirable to use universal multifractal model for anisotropic (i.e., different scaling in the vertical than in the horizontal direction) inhomogeneous turbulence (e.g., Lovejoy & Schertzer, 2010;Schertzer & Lovejoy, 2011). The anisotropic universal multifractal approach typically predicts the passive scalar turbulent H in horizontal to have values between 0.3 and 0.45 (De Montera et al., 2011;Lovejoy & Schertzer, 2010), consistent with SST universal multifractal spectra from Schmitt et al. (1996), , , Lovejoy et al. (2001), andDe Montera et al. (2010). This scaling is expected to be valid up to the scale of the largest eddies (∼150 km, Frenger et al., 2015). This range of H is not very far from the prediction of isotropic 3-D turbulence (H = 0.33) and is consistent with the prediction of surface quasi-geostrophic turbulence (Blumen, 1978). This means that the SST scaling profile with H in the 0.3-0.45 range is quite robust within multiple different models of turbulence, which makes it also harder to distinguish which of the models better represents reality. Above the largest eddy scale, there is no eddy-induced variability, so the scaling slope flattens (H = 0). We confirmed that the 0.3-0.45 turbulent pas- sive scalar range of H values is relevant for this study by calculating universal multifractal SST parameters for the EO data (see supporting information S1: Table S1) and using the approach of Schmitt et al. (1996) and De Montera et al. (2010Montera et al. ( , 2011. Using the EO SST data and the universal multifractal theory has advantage over analyzing model outputs, as we avoid some of the limitation imposed on eddies by the NEMO model resolution. A detailed handling of the universal multifractal methodology, together with the theoretical origin of power law scaling, can be found in section S2 of supporting information S1.  Note. As in equation (3), H 1 and H 2 are the two scaling exponents below (H 1 ) and above (H 2 ) the scale of the scaling break ( s ). The abbreviations for SST drivers are introduced in the section 1, and the reader is reminded of them in Figure 7 (explained in the caption). The table also indicates the scales of the scaling break (transition from H 1 scaling exponent to H 2 ). The scaling exponents H 1 and H 2 are rounded to two decimal places and s >100 km are rounded to multiples of 10 and s <100 km to multiples of 5. In some cases the parameters (i.e., s for ZWS and bathymetry) could not be on the 95% confidence level rounded to a unique value. In those cases we show the approximate 95% confidence interval. SST = sea surface temperature.

Results
We explored how accurately the different model resolutions reproduce the satellite EO spatial patterns of variability using bias-corrected root-mean-square error (RMSE; see supporting information S1, or, e.g., Taylor, 2001) as the skill metric ( Figure 4). First, Figure 4 shows that the finer-resolution models, at each scale, outperform the coarser resolutions (with the exception of the 1/4 • model in the Pacific). Second, and perhaps surprisingly, Figure 4 shows that the difference in how the finer-and the coarser-resolution models perform becomes more pronounced with increasing scale. This is because the finer-resolution model skill in representing EO SST spatial variability patterns improves more rapidly with scale (steeper slope in Figure 4) than the coarser resolution model skill. For example, Figure 4 shows that at the 75-to 750-km range of scales the 1/12 • resolution model improves in RMSE by 30-50%; the 1/4 • model by 20-30% and; and the 1 • model by <5%.
The numerical comparison between the model and the EO SST mosv has shown that each model resolution is at 4-5 times the model resolution scale in agreement with the observations within 10%. These results are shown in Figure 5, where the dashed lines represent 10% and 30% deviation from the EO SST mosv. Figure 5 shows that the model SST mosv always appears to be within the ±10% (dashed line) envelope around the EO SST mosv from a scale that is less or equal to the 4-5 times model resolution scale. At sufficiently large scales (O∼500 km) the different resolution models match the EO to within 3%.
EO SST mosv scaling can be described by two power law regimes (straight lines in log-log space) with a scaling break between 180 and 250 km (depending on the region, see Table 2 and Figures 6a, 6c, and 6e). Below 180-250 km the SST scaling slope is flatter, which means that variability reduces more slowly with the reduction in scale. The SST scaling slope becomes steeper above the 180-to 250-km scale, with H (corresponding to the scaling slope) close to 1 (see Table 2). H ≈ 1 means the SST variability scales above 180-250 km linearly with scale, i.e., that N% reduction in scale leads to N% reduction in SST variability (see equation (3), or equation (S4) in supporting information S1).
The scaling of SST drivers can also be split into two power law regimes separated by a scaling break s , with different s values for different SST drivers (Figures 7a, 7c, and 7e). For all the fields (SST and SST drivers) the piecewise power law fits the scaling of the mosv with an error <4% (in most cases it is <2%). Unlike SST, the SST drivers scale below s with steeper slopes than above s . Table 2 summarizes, for both SST and SST drivers, the scaling exponents and the scales for the scaling breaks. Table 2 demonstrates that (1) there is relatively little regional (Pacific, Indian, and Atlantic) difference in SST and SST drivers scaling.
(2) The scaling break occurs at a range of scales which can differ by an order of magnitude (30-800 km). We estimated the uncertainty of the parameter values presented in Table 2 from a suitably chosen random sample (of 1,000 members) generated by adding white random Gaussian noise to the mosv, with the mosv uncertainty used as standard deviation. Due to large number of used data (>10 6 ), the mosv uncertainty is at each scale only 0.1% of the mosv value (with 95% confidence). The uncertainty in the scaling exponents H 1 and H 2 was found to be very small (for each scaling exponent within the order of 10 −3 ). The uncertainty in Figure 8. (a, c, and e) The scaling of SST zonal, meridional, and total magnitude of spatial variability over 8-to 1,500-km range of scales. As in Figures 5-7 both axes are on a log scale. Y axis of the right-side panels (b, d, and f) shows the logarithm of meridional-to-zonal mean increment ratio for SST and SST drivers on the 8-to 1,500-km interval (log 10 ( ) shown in the x axis). The abbreviations for the SST drivers used in the right-side panels (plot legends) are explained in the caption to Figure 7. The horizontal lines in the right panels show the interval where the distributions have zonal increments within ±50% of the meridional increments, or equivalently the meridional-to-zonal increment ratio in the 0.67-1.5 range. SST = sea surface temperature.

10.1029/2018JC014791
s is in some specific cases nonnegligible (∼10 km), and in those cases Table 2 shows the approximate 95% confidence interval.
We analyzed the relative distribution of SST drivers mosv across the 8-to 1,500-km scales (Figures 7b, 7d,  and 7f). Based on Figure 7 one sees that resolving N = 85% of largest (1,500 km) scale SST driver subgrid variability will put the SST drivers resolution N at similar scales to the model resolution. We are therefore looking for the scale where the driver is left with 15% of its 1,500-km mosv. SST drivers can be split into three groups (Figures 7b, 7d, and 7f): (1) dominant large-scale variability drivers (SWR, ZWS), with N = 85% of variability resolved ( N ) at the 80-to 100-km scale; (2) medium-scale variability drivers (bathymetry and the AHF), with 85% variability resolved at the 20-to 40-km scales; and (3) small-scale variability drivers (WSC), with 85% of variability resolved between the 10-to 20-km scales. To resolve 85% of turbulent eddy induced 1,500-km SST subgrid variability, one needs an even smaller scale than 10 km. The SST turbulent eddy scaling exponents H were estimated from the universal multifractal exponents (see supporting information S1: Table S1) and the theory of Schmitt et al. (1996) and De Montera et al. (2010. They were found to be within the 0.33-0.4 interval, and this places the passive tracer turbulent SST N around the 1-km scale. There are a few points that merit a mention: (1) In the Indian region there is little difference between drivers from the second and third groups as both have resolutions N (N = 85%) near the boundary of the two intervals at the 20-km scale.
(2) Using equation (S8) from supporting information S1, one can evaluate the sensitivity of N to the choice of regional scale L = 1, 500 km. From equation (S8), the larger the absolute value of exponent the more sensitive N is to the uncertainty in L. Let us select AHF as reference driver (see section S1 of supporting information S1), the exponents can be calculated using the scaling exponents from Table 2 to be (on average for the three regions): 0.04 for ZWS (MWS), 0.16 for SWR, −0.29 for bathymetry, and −0.16 for WSC. A negative exponent means that N decreases if L increases and vice versa. Bathymetry has the largest absolute value of , which is approximately equal to the third root (≈0.33). This implies that decreasing the value of L by 30% would increase N for bathymetry by 12%. By comparison, the change to the N value corresponding to the 30% decrease of L is for the other SST drivers less than 6% (only 1.5% for zonal wind stress). Therefore, the classification of SST drivers resolution scales N is relatively unaffected by ambiguity in L.
(3) The scale of N ≈ 1 km is clearly indicative of the fact that turbulent eddies drive the small-scale SST variability. The exact value of N is however subject to some relatively large uncertainties: First, the model of Schmitt et al. (1996) and De Montera et al. (2010 is based on several assumptions, such as the closeness between the energy and tracer concentration variance flux multifractal parameters (De Montera et al., 2010). Second, in the SO eddies move with a velocity of roughly 22 km per week ; in a 5-day model output, they will as a consequence be smeared through a larger volume. To what extent this smearing affects the multifractal exponents (supporting information S1: section S2) is not known. Third, the ∼1-km variability scale is lower than the scale of EO SST data resolution, so the N value is based on extrapolating the EO SST power law scaling from the observed range to the scales of the order of ∼ 1 km. The turbulent N is also more sensitive to the ambiguity in L = 1, 500 km (turbulent exponent from supporting information S1, equation (S8), is = −1.1). However, none of these uncertainties is expected to change the order of magnitude difference in N between the turbulent passive tracer SST and the small scale resolution drivers (WSC).
SST spatial scaling is approximately isotropic (zonal-to-meridional increment ratio not far from 1) below the 50-km scale (Figure 8). At larger scales the zonal increments scaling flattens, while the meridional increment scaling becomes steeper. As shown in Figures 8b, 8d, and 8f the meridional increments eventually become 7-10 times larger than the zonal increments on the 1,500-km scale. Meridional increments also dominate over zonal increments for the (incoming) SWR and ZWS (Figures 8b, 8d, and 8f). For ZWS the meridional-to-zonal increment ratio changes only slightly with scale. However, for SWR the ratio doubles between the 8-and 1,500-km scale. For AHF the meridional increments are slightly larger (ratio 1-1.5) than the zonal increments, whereas bathymetry increments are (slightly) more zonal (ratio 0.67-1). The uncertainty for the meridional-to-zonal ratio was estimated using the method applied to the piecewise power law parameters and it was found to be very small (in all cases < 0.3% from the ratio).  (2)) versus logarithm of time for Pacific (a), Indian (c), and Atlantic (d). The different lines are the spatially low-pass-filtered EO SST data at a range of scales varying between 20 and 1,600 km. The 8-km data are the unfiltered EO SST data set (8 km is the EO SST spatial resolution scale). The dashed curves mark 98% and 95% of the unfiltered data motv and the vertical dashed lines the time scales where the motv of the filtered data matches the motv of the unfiltered data within 2% (the colors of the vertical lines match with the color of the spatial low-pass filtering). (b) Pacific log-time (log( )) versus log scale (log( )) relationship between the spatial scale of low-pass filtering ( ) and the time scale ( ) when the processes filtered from the data account for less than 2%, 3%, and 5% of motv. The purpose of panel (b) is to show the connection between the process spatial and temporal scales found in the SST data. The vertical dashed line shows the scale of the spatial scaling break (180 km) found in the EO SST magnitude of spatial variability scaling in the Pacific. Similar ( ) curves can be reproduced for the other two regions (Indian and Atlantic). To avoid the figure becoming crowded, only the most interesting, Pacific, case is shown. SST = sea surface temperature; EO = Earth Observation.
Finally, we analyzed the EO SST motv (see equation (2)) for different spatially low-pass-filtered data (filtered on the 20-to 1,600-km range of scales; Figure 9). The spatial low-pass filter at the scale removes all the processes that have impact on SST below . At the time scale ( ) where these removed processes are no longer important, the SST motv (equation (2)) of the spatially low-pass-filtered data becomes roughly equal to the motv of the unfiltered (original) data with the ∼8-km resolution. By defining a threshold for when the time variability of the filtered data is sufficiently close to the time variability of the unfiltered data, we can calculate for each length scale its associated time scale . The spatial ( ) and temporal ( ) process scales have been matched (Figure 9b) using thresholds equal to 2%, 3%, and 5% of the unfiltered data motv. Although the spatio-temporal relationship ( ) has a certain ambiguity (due to the ambiguity in the value of the threshold), Figure 9b shows that the scale of the SST spatial scaling break ∼ 200 km (see Figure 6) corresponds to a scale where the ( ) curves flatten. The implications of this will be discussed later. Figures 9a, 9c, and 9d also show that at the 1,600-km spatial filter scale (Pacific) and at ≥800-km scale (Indian and Atlantic) motv fails to converge (within the half-year time period) to the unfiltered data (with a gap of at least 5% between the two curves). This is because the low-pass filtering removed at ∼1,000 km also supraseasonal processes, including the long-term SST SO geography (e.g., 10 5 -year-averaged spatial distribution).

Summary and Discussion
Each of the questions raised in section 1 will now be taken in turn for summary and discussion:

Question 1: How does model resolution impact on model skill in representing SO SST above the submesoscale?
In summary, 1. increasing model resolution has a substantially positive impact on how the model represents EO spatial variability patterns, and this impact increases with scale ( Figure 4). 2. At the model resolution scale the model magnitude of SST spatial variability is around 15-30% lower than the same scale EO mosv (the difference is typically close to 30%; see Figure 5). 3. From 4-5 times the resolution scale all three model configurations represent the magnitude of EO spatial variability within ±10% (Figure 5). 4. The accuracy of how the model represents the EO mosv improves with scale, and above 500 km each of the models matches the EO to within ±3% ( Figure 5).
The model has limited capability to correctly represent advection of tracers (SST) near its grid scale, since it has to parametrize subgrid eddy flux by the grid-scale viscosity/diffusion coefficient. At scales sufficiently larger than the model resolution the model will represent SST better, as at those scales the grid-scale viscosity/diffusion no longer impact on the advection of tracers (e.g., Su et al., 2016aSu et al., , 2016bSu et al., , 2018. Figure 4 indeed shows that the RMSE skill improves between the model resolution scale and the 750-km scale by 5-55%, depending on the model resolution. As expected, the finer-resolution model RMSE skill is always better than the coarser resolution model (Figure 4, with the exception of 1/4 • model at its resolution scale in the Indian region). However, we also expected that the difference in skill between different model resolutions should become smaller with scale, since the impact of model resolution will gradually diminish with scale. Counterintuitively, Figure 4 shows the opposite: the difference in skill (RMSE) between the finer-resolution model and the coarser-resolution model grows with scale. While at the 1 • model resolution scale (≈70 km) the 1/12 • model outperforms (reduces RMSE) the 1/4 • model by ≈12% and the 1 • model by ≈30%, at the 750-km scale (∼12 • ) the 1/12 • model outperforms the 1/4 • model by ≈30% and the 1 • model by ≈60%.
Using magnitude of SST spatial variability as skill, it can be concluded ( Figure 5) that all three model resolutions represent the EO SST with sufficient accuracy (within 10%) from 4-5 model resolution scale. This is consistent with the recent observation from Soufflet et al. (2016). While the magnitude of EO SST variability has very similar scaling in all three regions (Table 2), the model SST scaling (especially the 1/12 • model) has important differences between regions ( Figure 5). In the Pacific the 1/12 • model scales between 50 and 200 km similarly to the EO, whereas for both the Atlantic and Indian regions, the model has on the 50-to 200-km scales flatter scaling than the EO data. In these same two regions the model overestimates on the 50-to 200-km scales the magnitude of SST variability (although still within ±10%). In the Atlantic and Indian regions, the 1/12 • model matches EO within ±3% from approximately 200-km scale (30 times model resolution), whereas in the Pacific it matches EO with the same accuracy at a 50-km scale ( Figure 5).
It is difficult to explain the regional differences in the 1/12 • model SST mosv scaling below the 180-km scale. The larger (than EO) 1/12 • model 30-to 80-km scale mosv in the Atlantic and Indian regions could indicate that OI has reduced the EO mosv below ≈100-km scale (the OI smoothing radius). However, the relative effect of OI smoothing on the scaling of EO mosv is expected to gradually increase with the decreasing scale. This would appear as a sudden reduction in SST mosv close to the EO SST resolution scale, in a similar fashion to what happens to the 1/4 • and 1 • models at their resolution scales ( Figure 5). Contrary to this expectation the magnitude of EO mosv is scale invariant below the 100-km scale. Moreover, the scale invariance of EO SST mosv on ∼8-100 km is consistent with other scaling analyses of IR EO SST products from the North-West European Shelf, which included higher spatial and time resolutions (Skákala et al., 2016). Even the SST scaling exponents (H 1 ) found here on the 8-to 100-km interval are similar to the values from Skákala et al., 2016Skákala et al., (2016 H 1 around 0.6-0.8). These are strong arguments to suggest that spatial smoothing introduced by time averaging, or OI had a relatively small effect on small scale SST mosv. Figure 5 shows that the 1/4 • and 1 • models have, at their resolution scale, about 30% lower magnitude of SST variability than the EO (with the exception of 1/4 • model in the Atlantic). There are two reasons for this. First, the mosv of the finer-resolution EO data set at the 1/4 • (or 1 • scale) also includes intrapixel variability. This is because among two arbitrary 1/4 • (1 • ) model pixels separated by scale the 1/12 • model connects multiple 1/12 • pixels (separated by the same scale ). Second, models effectively average out the dynamics close to their resolution scale and this also reduces the SST mosv. We determined the magnitude of missing intrapixel variability by averaging out all the intrapixel variability of the finer-resolution model on the coarser-resolution model scale. The SST mosv (equation (1)) was then compared on the coarser resolution scale with the finer-resolution model. Using this procedure, we calculated that the missing intrapixel variability lowers the coarser resolution model grid-scale variability by ∼10% and therefore contributes slightly less than half to the missing SST mosv near the model grid scale.
Question 2: Which model input drivers need to be well resolved in order to reproduce the SO SST variability at a given scale?
In summary, 1. the optimal resolution N is defined as the scale resolving N=85% of the SST driver 1,500-km subgrid variability. 2. Based on the optimal resolution scale N , we classified the SST drivers (see Figure 7, also supporting information S1: Figure S7) as follows: 2.1. dominant large-scale variability drivers (short-wave radiation (SWR), wind stress vector components, and cubic wind speed), with a resolution scale ( N ) between 80-and 100-km scales; 2.2. dominant medium-scale variability drivers (bathymetry and atmospheric heat fluxes (AHFs)) with N on the 20-to 40-km scales; 2.3. dominant small-scale variability drivers (wind stress curl (WSC)), with N between 10 and 20 km; and 2.4. in order to capture 85% of eddy-induced (1,500 km) SST variability a model resolution below 10-km scale is required (≈1 km).
A great deal of this is not surprising. SWR and ZWS (MWS) are responsible for the seasonal forcing and the SST geography. This implies SWR and wind stress drive variability on the regional (∼1,000 km) spatial scales and are therefore expected to have a dominant large-scale variability. The AHFs and bathymetry (medium-scale SST input drivers) affect all the analyzed scales, with the smaller-scale AHF variability mostly due to the fast atmospheric response to SST perturbations and the fast atmospheric dynamics (Hausmann, Czaja, et al., 2016). The small variability scales for turbulent passive scalars are linked to the limited SO eddy size (≲150 km; Frenger et al., 2015). Similarly, it is well known (e.g., Frankignoul, 1985;Haney et al., 1978;Tamsitt et al., 2016) that WSC has almost no impact on SST on synoptic scales, but potentially large impact on the small-scale SST variability. For example, to represent adequately the Benguela upwelling system, one needs high-resolution atmospheric model and eddy resolving ocean model (Small et al., 2015). Although the classification of the SST input driver-dominant variability scales can be anticipated from Figure 1, the scaling analysis enables us to precisely quantify the relationship between the different model SST input driver optimal resolutions.

Question 3: What does the scaling of SO SST variability tell us about SST dynamics above the submesoscale?
In summary, 1. the SST scales with two power law regimes (equation (3)) separated by a scaling break s at 180-to 250-km scale (depending on the region; see Figure 6 and Table 2). 2. The SST scaling break (∼200 km) is a scale of separation between two SST characteristic dynamical regimes captured by the two distinct power laws: (1) the subseasonal processes (mesoscale eddies and fronts, Ekman pumping, latent, and sensible heat fluxes) and (2) the seasonal and supraseasonal modes that depend mainly on the short-wave radiation and (zonal) wind stress. (This is supported by the analysis presented in Figures 1, 8, and 9.) It has been shown in Figure 6 that SST scales as a piecewise power law with a scaling break ( s ) at ∼200-km scale. This result can be interpreted using the diagram in Figure 1, which suggests that s is a separation between two characteristic dynamical regimes: the subseasonal SST dynamics associated with eddies, fronts, air-sea heat flux, and Ekman pumping, and the seasonal and supraseasonal dynamics associated with planetary scale forcing. This interpretation of the SST scaling law is further supported by Figure 9, showing the connection between spatial and temporal scales of the observed SST processes. Although the precise connection between the spatial ( ) and temporal ( ) scales is to some extent ambiguous, as it depends on the threshold below which a process is neglected, the shape of the ( ) curve is independent of the threshold value ( Figure 9). The SST scaling break s corresponds to the spatial scale when the ( ) curves flatten. A steep ( ) curve means that spatial filtering has a large impact on the SST time variability curve, because it removed substantial number of processes associated with the SST dynamics. A flat ( ) curve implies that no additional processes have been removed by the spatial filter. s then corresponds to the spatial scale where the numerous subseasonal SST processes vanish and for the next range of spatial scales the SST is uniformly governed by the seasonal forcing. However, Figures 9a, 9c, and 9d also show that at sufficiently large spatial scales (∼1,500 km in Pacific and ∼800 km in Indian and Atlantic) the SST variability is dominated by supraseasonal modes. The supraseasonal modes can be anything from interannual climate modes (see Figure 1) to the SST geography (e.g., 10 5 -year time average for each spatial location). Analysis presented in supporting information S1: Figure S8 indicates that from the supraseasonal modes the interannual climate modes such as El Niño-Southern Oscillation (ENSO) or Southern Annular Mode (Kostov et al., 2017;Verdy et al., 2006;Yang et al., 2007) have a relatively secondary effect on the scaling of SST mosv. Comparison of the 2010 and 2012 EO SST data (see supporting information S1: Figure S8) shows the maximum difference in the SST mosv on the level of 1%, which is barely statistically significant. This is negligible compared to the large reduction of SST variability (on the level of tens of percent; see Figure 9) when the spatial filter was applied to ∼1,000-km scales. It is very likely that the supraseasonal modes observed in Figure 9 at ∼1,000 km scales can be primarily attributed to SST geography. The major contributors to SST geography are air-sea fluxes (both SWR and AHFs), Ekman transport, and the heat advection by geostrophic currents (see the 5-year analysis of SO heat budget in (Tamsitt et al., 2016(Tamsitt et al., ), 2016. Interestingly, the supraseasonal modes influence SST spatial variability on smaller scales in the Indian and Atlantic regions than in the Pacific (Figure 9). This is consistent with the analysis of Tamsitt et al. (2016), suggesting that this smaller-scale SST variability is due to heat transport by geostrophic eddies and vertical transport. Curiously, part of this small-scale SST variability could be also attributed to the observed larger climate variability in the Indian and Atlantic regions than in the Pacific (see Deser et al., 2010).
We can analyze the two SST scaling regimes also from a slightly different perspective. Let us start with a simplifying assumption that the impact of SST drivers on SST variability depends mainly on the magnitude of SST driver variability. We then ask which SST drivers have (among all the SST drivers considered) the largest percent of 1,500-km scale variability at the scales above 300 km (range of scales above the scale of the scaling break and well separated from it). Figure 7 shows that these are wind fields (ZWS and MWS) and SWR with ≈60% of variability above the 300-km scale. These are the SST drivers largely responsible for seasonal and supraseasonal dynamics (including SST geography). For comparison the turbulent passive tracer SST has no additional variability above 300 km, WSC and bathymetry have around 30%, and AHF has around 40%. Conversely below the 100-km scale (below the scales of scaling break and well separated from it), passive tracer (turbulent) SST has 86% from its 1,500-km scale variability, WSC around 45%, bathymetry and AHF around 25-30%, while SWR and ZWS only 15% (Figure 7). This suggests that wind fields and SWR are the most dominant drivers for SST variability at the scales above 300 km, whereas turbulent eddies are the most dominant SST driver below 100-km scale. The large impact of SO turbulent eddies on SST signature on scales up to 150 km is well known in the literature (Chelton, 2013;Foppert et al., 2017;Frenger et al., 2015Frenger et al., , 2013Mazloff et al., 2010;Olbers et al., 2004;Rintoul & Sokolov, 2001). Even in the southern parts (south of 50 • S) of the Indian and Pacific SO regions the eddies reach scales of several tens of kilometers Frenger et al., 2015;Hausmann, McGillicuddy, et al., 2016, which is captured by the EO ≈1/12 • resolution. The difference in the anticyclonic and cyclonic SST anomaly in the SO can be >1 • C . It is therefore a plausible assumption that the relatively flat scaling slope of the turbulent SST (0.33-0.41) is responsible for the flatter SST scaling slope at scales below s (the scale of the SST scaling break), when compared to the SST scaling slope above s . Above the scaling break, the SST scaling steepens, consistently with SWR and ZWS scaling; however, curiously, when it comes to the SST scaling close to the maximum L = 1, 500-km scale, SST has a steeper scaling slope than SWR and the wind fields. It would be interesting to further explore whether this results from the nonlinearity of the SST driver impact on the SST variability as suggested by the analysis in supporting information S1 (see the equations (S18-S21) and therein).
The observed scales of model input driver impact on SST are also seen in the analysis of SST and SST drivers meridional-to-zonal increment ratio (Figure 8). Figure 8 indicates that there are two SST drivers with a large degree of horizontal anisotropy: ZWS (over the whole 8-to 1,500-km range of scales) and SWR (mostly above the ≈200-km scale), both of those SST drivers having meridional increments much larger (by more than 50%) than the zonal increments. The SST meridional-to-zonal increment ratio from Figure 8 suggests that while 10.1029/2018JC014791 SST distribution is isotropic on ∼10-km scales, it becomes quickly anisotropic on ∼100-km scales ( Figure 8) and the meridional SST increments at the 1,500-km scale become a factor of 7-10 times larger than the zonal increments. Wherever SST distributions are horizontally isotropic (∼10-km scales), we assume that the distribution is shaped by horizontally isotropic SST drivers (turbulence, to some degree also bathymetry, WSC, AHF, but at ∼10-km scales it can be also SWR). Wherever SST distributions are strongly anisotropic (≿100 km), it will become dominated by strongly anisotropic SST drivers (ZWS, SWR). The SST transition from isotropy to anisotropy does not seem to be driven by a specific SST driver increase in anisotropy (see Figure 8, SWR has a substantial increase in anisotropy from larger, 200-km scales). The SST transition from isotropy to anisotropy is more likely to be a consequence of strongly anisotropic SST drivers (ZWS and SWR), responsible for the seasonal and supraseasonal SST dynamics, becoming at the ∼100-km scales dominant over SST isotropic drivers in their impact on SST distribution.
Although this paper has studied how SST drivers influence SST scaling, within the marine boundary layer there is substantial feedback of SST on atmospheric SST drivers (Chelton, 2013;Chelton & Wentz, 2005;Chelton & Xie, 2010;Chelton et al., 2004;Frenger et al., 2013;Xie, 2004). For example, Frenger et al. (2013) have found evidence that there is a significant link between eddy induced SST variability and surface winds at the mesoscale: colder SST associated with cyclonic mesoscale eddies weakens surface winds, while the winds are strengthened by warm SST associated with anticyclonic eddies (see also discussion in Chelton, 2013). SST also has significant impact on low-level clouds (influencing incoming short-wave radiation) and precipitation (Chelton, 2013;Frenger et al., 2013). It is then possible that in reverse one can associate some SST driver scaling break (see Table 2) with a maximum scale of SST impact on atmospheric fields in the marine boundary layer (Chelton, 2013;Frenger et al., 2013).

Conclusions
In this paper we have introduced a methodology based on scaling analysis techniques to (1) evaluate model resolution skill in representing SO SST, (2) evaluate the resolution needed to capture the variability of model input SST drivers, and (3) learn about dominant processes responsible for SST dynamics across a wide range of spatial scales (above the submesoscale).
We have found that the scale s ∼200 km provides important separation between two distinct dynamical regimes: below s SST is governed by the subseasonal dynamics (eddies, fronts, air-sea heat flux, and Ekman pumping), while above s the SST is driven by seasonal forcing and longer modes. The connection between temporal and spatial scales is particularly important: it could help inform, for example, filtering steps when wanting to separate eddy and mean (i.e., seasonal and supraseasonal) signals for cases when only time or spatial content is available. Here we have in mind a separation time scale useful for Reynolds decomposition. Quantifying the spatial separation between the subseasonal and longer-term processes is also important from the modeling perspective, as it informs what processes need to be well represented by the model across a wide range of spatial resolutions.
This paper necessarily concentrated on multiple analyses of a single year (2010) due to data volume and availability issues. The scaling analysis applied to 2012 data showed little impact of intra-annual modes (ENSO and Southern Annular Mode) on the conclusions presented in this study. However, better exploration of what happens on the climatic time scales will be desirable. For example, it would be interesting to determine whether the scaling of SST spatial variability is influenced by the climate change and its impact on the SO eddy kinetic energy (Hogg et al., 2015;Meredith, 2016;Meredith & Hogg, 2006).