Global Impacts of Subseasonal (<60 Day) Wind Variability on Ocean Surface Stress, Buoyancy Flux, and Mixed Layer Depth

Subseasonal surface wind variability strongly impacts the annual mean and subseasonal turbulent atmospheric surface fluxes. However, the impacts of subseasonal wind variability on the ocean are not fully understood. Here, we quantify the sensitivity of the ocean surface stress (τ), buoyancy flux (B), and mixed layer depth (MLD) to subseasonal wind variability in both a one-dimensional (1-D) vertical column model and a three-dimensional (3-D) global mesoscale-resolving ocean/sea ice model. The winds are smoothed by time filtering the pseudo-stresses, so the mean stress is approximately unchanged, and some important surface flux feedbacks are retained. The 1-D results quantify the sensitivities to wind variability at different time scales from 120 days to 1 day at a few sites. The 3-D results quantify the sensitivities to wind variability shorter than 60 days at all locations, and comparisons between 1-D and 3-D results highlight the importance of 3-D ocean dynamics. Globally, subseasonal winds explain virtually all of subseasonal τ variance, about half of subseasonal B variance but only a quarter of subseasonal MLD variance. Subseasonal winds also explain about a fifth of the annual mean MLD and a similar and spatially correlated fraction of the mean friction velocity, u∗ = √|τ|∕ρsw where ρsw is the density of seawater. Hence, the subseasonal MLD variance is relatively insensitive to subseasonal winds despite their strong impact on local B and τ variability, but the mean MLD is relatively sensitive to subseasonal winds to the extent that they modify the mean u*, and both of these sensitivities are modified by 3-D ocean dynamics.


Introduction
Several crucial ocean and Earth system processes, including water mass formation and circulation (Hanawa and Talley, 2001;Stommel, 1979), air-sea exchange (Frankignoul and Hasselmann, 1977;Kraus and Turner, 1967;Rodgers et al., 2014;Stevenson and Niiler, 1983), and biogeochemistry (Fasham et al., 1990;Fauchereau et al., 2011;Sverdrup, 1953), are highly sensitive to the ocean surface mixed layer depth (MLD). Due to the rapid increase in Argo observations over the last decade, the global seasonal climatology (including mean and variance) of the MLD is increasingly well observed on the global scale (Holte et al., 2017). In addition, it is well known that the ocean surface stress and buoyancy flux inject kinetic energy or buoyancy/stratification to form, maintain, or change the MLD (Deardorff et al., 1969;Kraus and Turner, 1967;Kato and Phillips, 1969;Large et al., 1994;Pollard et al., 1972;Price et al., 1986;McWilliams et al., 1997). However, the dynamics of the climatological MLD are not fully understood and remain challenging to model globally (Belcher et al., 2012;D'Asaro, 2014;DuVivier et al., 2018).
Models can fill in observational gaps by providing more comprehensive spatiotemporal coverage and deeper insights into the processes controlling the MLD. Previously, mostly one-dimensional (1-D) vertical column models have been used to build understanding of the ocean mixed layer's response to a synoptic wind event (Pollard et al., 1972;Niiler, 1975), the diurnal cycle (Price et al., 1986), the seasonal cycle (Gill and Turner, 1976;Martin, 1985;Thompson, 1976), and the longer-term equilibrium response to forcing Price, 1998;Price and Sundermeyer, 1999). Turbulence-resolving process models have been used to validate 1-D model mixing parameterizations on time scales of a week or less (Large and Gent, 1999;Van Roekel et al., 2018), and observations have been used to verify basic properties of model solutions on longer time scales. Such studies, in combination with observations, have shown that interannual variations in subseasonal atmospheric variability or even just a few crucial synoptic events can substantially influence the seasonal evolution of the MLD at some subpolar sites (Large et al., 1986;Waniek, 2003), while subseasonal atmospheric variability may have a far smaller, though still significant, impact on the MLD at lower latitudes (Bernie et al., 2005;Shinoda and Hendon, 1998;Zhou et al., 2018).
No previous study has presented a global assessment of the impacts of subseasonal wind variability on the MLD using a 1-D model, because such models can be inadequate in some regions. For example, lateral advection by mean currents, which must be parameterized in a 1-D model, can contribute significantly to setting the upper-ocean temperature, hence heat flux, stratification, and the mean MLD in some regions, such as western boundary currents (Niiler, 1975). In addition, mesoscale and submesoscale oceanic variability, which must also be parameterized in a 1-D model, can significantly impact the mean and seasonal cycle of the MLD in some regions, like the Subantarctic Zone of the Southern Ocean (du Plessis et al., 2017;DuVivier et al., 2018;Fox-Kemper et al., 2011;Lee et al., 2011;Li and Lee, 2017;Swart et al., 2015). Indeed, high-resolution simulations that resolve some subseasonal atmospheric variability and ocean mesoscale variability show that (among other properties; (e.g., Delworth et al., 2012;Griffies et al., 2015;Kirtman et al., 2012;Maltrud and McClean, 2005;Maltrud et al., 2008Maltrud et al., , 2010McClean et al., 2011;Small et al., 2014;Winton et al., 2014) the mean MLD depends on the horizontal grid spacing of the atmosphere and ocean model at resolutions between 2 • and 0.1 Harrison et al., 2018;Lee et al., 2011;Sein et al., 2018), which suggests that both atmospheric and oceanic mesoscale processes (which vary on subseasonal time scales) have a substantial impact on the annual mean MLD. Hence, high-resolution three-dimensional (3-D) ocean circulation models are more suitable than 1-D models for simulating the dynamics of the climatological mean, seasonal cycle, and subseasonal variability of the MLD on a global scale.
Here, we quantify the sensitivity of the climatological annual mean, seasonal cycle, and subseasonal variance of the MLD to subseasonal winds, via their effects on ocean surface stress and buoyancy flux, at the finest resolved spatial scales (3-11 km) in a global mesoscale-resolving ocean/sea ice numerical circulation model for the first time. Although it is well known that subseasonal winds have a significant impact on the stress, buoyancy flux, and MLD, no previous study has quantified these impacts in a global model that resolves mesoscales, which generally dominate upper-ocean variance (Delworth et al., 2012;Ferrari and Wunsch, 2009;Kirtman et al., 2012;Small et al., 2014;Wunsch, 2002) and have a significant impact on the mean and seasonal cycle of the MLD (DuVivier et al., 2018;Hausmann et al., 2017;Harrison et al., 2018;Lee et al., 2011;Sein et al., 2018). In this modeling context, the motivating questions for this study are as follows: What fraction of the climatological mean, seasonal cycle, and subseasonal variability of the MLD is attributable to subseasonal winds, and how do these fractions vary across the globe? What fraction of MLD variance is subseasonal, and how does subseasonal MLD variance vary across the globe? Moreover, to what degree does 1-D vertical mixing physics and local surface fluxes or 3-D ocean dynamics control the sensitivity of the MLD to subseasonal winds? Can the sensitivity of the MLD to subseasonal winds be explained by changes in the local surface fluxes of buoyancy and momentum using theoretical scaling arguments?

Model Setup 2.1.1. 3-D Ocean/Sea Ice Model
To assess the role of subseasonal surface wind variability in setting climatological mean, seasonal cycle, and subseasonal variability of the MLD, we present two 5-year simulations of global ocean and sea ice dynamics in a nominal 0.1 • resolution configuration of the Community Earth System Model (CESM), which is forced by atmospheric surface conditions and runoff that are derived from data sets as described in section 2.2. The control (CTL) and "low-pass" (LP) simulations differ only in the specification of the surface wind. In LP, the wind is smoothed in time to remove variability on time scales shorter than 60 days as described in section 2.2. Differences in a given output field between simulations (CTL minus LP) quantify the effects of the subseasonal winds on the output in question.
The ocean model is the Parallel Ocean Program Version 2 (Smith et al., 2010) with 62 vertical levels ranging in thickness from 10 m over the top 150 m to 250 m at the bottom, and the sea ice model is the Community Ice Code Version 5 (Bailey et al., 2018). The equations are solved on a tripolar orthogonal curvilinear grid with poles in North America, Asia, and Antarctica. The horizontal grid contains 2,400 by 3,600 points that are spaced with an approximate cos(latitude) dependence thereby achieving a resolution that ranges from about 11 km at the equator to 3 km in the polar oceans. Hence, the grid resolution is sufficiently fine to resolve the first baroclinic Rossby radius of deformation, which ranges from about 100 km at 10 • latitude to 10 km at 60 • latitude (e.g., Chelton et al., 1998), over most of the global open ocean. Therefore, the model configuration is called "mesoscale resolving" (Hallberg, 2013). The configuration of the ocean model is similar to recent high-resolution ocean/sea ice configurations of CESM in most other ways as well (e.g., Bryan and Bachman, 2015;Delman et al., 2018;DuVivier et al., 2018;Harrison et al., 2018;Johnson et al., 2016;Soares et al., 2019). For example, as in these previous versions, runoff is mapped to a broad swath of grid cells near the coast, there are no explicit overflow, tidal mixing, near-inertial mixing, submesoscale, or mesoscale parameterizations, time stepping is achieved using the time-averaged leap frog scheme with a time step of 3.6 min, subgrid-scale turbulent vertical transport of tracers and momentum is computed with the K-profile parameterization (KPP) of Large et al. (1994), and subgrid-scale lateral transport is parameterized via a spatially variable biharmonic viscosity and diffusivity that scale with grid size as in Bryan and Bachman (2015). The ocean coupling interval is shorter than some previous runs with CESM, 30 min versus 6 or 24 hr, which is important for appropriately resolving inertial frequency ice-ocean responses (Roberts et al., 2015) (see also appendix C.6 of Griffies et al., 2009).
The initial conditions for the simulations presented here are from an ocean/sea ice model state in a previous 20-year run that used an earlier version of the CESM model code (documented in Bryan and Bachman, 2015). This previous run originally started from rest and an observational climatology of temperature and salinity, and it was executed on the same 0.1 • tri-pole grid with the same atmospheric forcing as in the CTL simulation (section 2.2.1).
To allow the solution to adjust to minor changes in the new model code/configuration, a spin-up was run for 22 months under CTL forcing from 1 February during Year 16 of the previous run by Bryan and Bachman (2015). On 1 December of Year 17, two branches are made: One continues under CTL forcing, and one is forced thereafter by LP forcing (section 2.2.2). The first month is excluded from the analysis, and Years 18-22 (inclusive) of these branches are analyzed and compared in section 3.
Experience shows that the upper-ocean drift decays rapidly over the first couple of decades of these simulations. Hence, the 18-year spin-up is important to help reduce (but not eliminate) model drift and isolate the signal associated with the change in forcing. Some measures of the drift are included in the supporting information. These results show that the 5-year drift in the simulated global mean surface kinetic energy, MLDs, temperature, and salinity is small relative to the amplitude of the seasonal cycle or the magnitude of the difference between the CTL and LP simulations ( Figure S4). Hence, we treat the five simulated years as essentially independent realizations without a trend in the analysis below, despite the fact that the model state is not perfectly in equilibrium.

1-D Ocean Model
One hypothesis is that changes in local atmospheric surface conditions and 1-D vertical mixing physics are responsible for a significant fraction of the sensitivity of the stress, buoyancy flux, and MLD to subseasonal winds. Here, we investigate this hypothesis and interpret the 3-D results with the aid of a series of 1-D simulations of the ocean surface fluxes and vertical mixing physics, without advective tendencies and pressure gradient forces, at a few example sites. In addition, the 1-D simulations complement the global spatial perspective of the 3-D simulations by showing how the various different time scales of subseasonal wind variability, ranging from 120 days to 1 day, impact the stress, buoyancy flux, and the MLD.
The 1-D model is implemented on a horizontally periodic 2-by-2 point grid in the Regional Ocean Modeling System (ROMS) (Haidvogel et al., 2008;Shchepetkin & McWilliams, 2005), which uses the KPP mixing scheme as implemented in the public repository myroms.org with KPP parameters modified from the ROMS defaults to match CESM. In addition, ROMS is modified to use the same bulk flux routine to calculate fluxes as a function of the atmospheric and oceanic states as in CESM (Large and Yeager, 2004;Large and Yeager, 2009;Small et al., 2015). For the sake of simplicity, Jerlov Water Type II is used to describe the penetration of solar radiation in all the 1-D simulations (Paulson and Simpson, 1977). The horizontal distance between each grid point is 10 km, and the time step is set to 20 min. The vertical grid includes 150 levels that are stretched over the full ocean depth using the stretching functions defined by equations (2.2) and (2.4) of Shchepetkin and McWilliams (2009) with h c = 250 m, s = 7, and b = 2 (e.g., the grid spacing ranges from Each 1-D simulation is of a deep, open ocean site defined by a box with size of 3 • by 5 • . The atmospheric forcing is defined by interpolating the atmospheric fields to a 40 km square patch in the middle of the box to avoid suppressing atmospheric variability via spatial averaging. Otherwise, the 1-D model is forced by the same atmospheric fields as the 3-D model, which are described in section 2.2.2, and the winds are modified in perturbation experiments using the same smoothing procedure as in the 3-D LP experiment, which is described in section 2.3.2. However, several different smoothing time scales, ranging from 120 days to 1 day, are applied to the winds. In one additional experiment, the wind variability with periods shorter than 5 days is amplified by multiplying the pseudo-stresses at those high frequencies by a factor of 3 (called 3xHF). The 1-D model is initialized in mid-February using the vertical profiles of temperature and salinity from the monthly climatology of the 3-D CTL simulation (horizontally averaged over the box), and the vertical profiles of temperature and salinity are restored over a time scale of 1 year toward that monthly climatology throughout each 1-D simulation. Velocity profiles are restored toward zero with a time scale of 10 days to prevent the accumulation of barotropic inertial energy, which is resonantly excited by wind stress variance at the local Coriolis frequency f but can only be damped by bottom friction in a 1-D model without momentum restoring. The restoring time scale of 1 year for temperature and salinity was chosen via test simulations with different restoring time scales at the Labrador Sea deep convection site (not shown), where it was found that some restoring was necessary to obtain a stable and reasonable repeating annual cycle in the 7-year integrations. The restoring time scale of 10 days for velocities was chosen to be reasonably representative of the observed decay time scale of inertial oscillations in the mixed layer (Park et al., 2009). The fifth full annual cycle is presented below, but the solutions have approximately converged so that the MLD is essentially a repeating annual cycle by this time.

Forcing 2.2.1. CTL Experiment
The surface atmospheric state in the CTL simulation is derived from the repeating "normal year" of the Coordinated Ocean Ice Reference Experiments (CORE) (Griffies et al., 2009;Large and Yeager, 2004), which consists of a single statistically typical annual cycle representative of the Years 1958-2000 (see section 2.3.1; section 9 of Large and Yeager (2004) as well as Griffies et al. (2009)). In the CTL simulation, the default CORE normal year protocol is followed. As prescribed in the default CORE protocol, sea surface salinity is restored toward observations based on the 1998 World Ocean Atlas with revisions in the Arctic (see section 4.7 of Large and Yeager (2004) and section B.3 of Griffies et al. (2009)). Certain aspects of the solutions (e.g., the meridional overturning circulation) are expected to be sensitive to this salinity restoring time scale in longer simulations (e.g., 500 years, as in Danabasoglu et al., 2014;Griffies et al., 2009), and the "optimal" restoring time scale for our simulations is unknown. Restoring time scales of 1-4 years have typically been used with CESM; here we use a 1-year restoring time scale in both simulations.
The ocean surface MLD is modified by surface fluxes of heat, freshwater, and momentum, which are calculated from the atmospheric, oceanic, and sea ice states using bulk formulas (Large and Yeager, 2004;Large and Yeager, 2009). For example, the stress on the ocean at the air-sea interface (and analogously at the ice-ocean interface) is calculated as where u is the ocean surface velocity, |u a | is the 10-m wind speed, u a = (u a , v a ) is the 10-m (zonal, meridional) wind vector, a is the air density, and c d is a variable exchange coefficient that is solved for. As described in Large and Yeager (2004), the atmospheric state variables are derived primarily from the NCEP reanalysis (Kalnay et al., 1996). Additionally, the wind speed and direction derived from the reanalysis are corrected by applying spatially variable correction factors, which do not vary with time, during the model execution (as described in Large and Yeager, 2004). The daily components of the ocean surface stress = ( x , y ) are saved during both the CTL and LP runs, and the magnitudes of these daily means | | = √ 2 x + 2 are calculated offline and analyzed in section 3. Here, and throughout the paper, the subscripts x and y indicate zonal and meridional directions, respectively.
The model employs qualitatively similar wind-speed-dependent expressions for the latent and sensible turbulent heat fluxes as well as evaporative turbulent freshwater flux (Large and Yeager, 2004;Large and Yeager, 10.1029/2019JC015166 2009. For example, the turbulent sensible heat flux is calculated as where a is the atmospheric potential temperature at 10 m, c p is the specific heat of seawater, c h is a variable exchange coefficient that is solved for, and SST is the sea surface temperature. All the turbulent fluxes come together with fluxes due to radiation, sea ice tendencies, precipitation, and runoff, to impact the MLD via the surface buoyancy flux Here, the freshwater contribution to the buoyancy flux is expressed as a virtual salt flux in the ocean model, and a positive buoyancy flux means that the ocean gains buoyancy and the surface water becomes less dense, either due to net freshwater input  surf or net heat input  surf . As described in Large and Yeager (2004) and Griffies et al. (2009), the atmospheric forcing fields, such as wind speed, humidity, and temperature, have a temporal resolution of 6 hr and a spatial resolution of about 2 • . Shortwave and longwave radiative forcing are provided as daily averages. Precipitation forcing is provided as monthly averages. These are all linearly interpolated to the model time steps internally via the coupling infrastructure. The thermal expansion coefficient (which ranges from −0.3 × 10 −4 to 3.4 × 10 −4 K −1 ) and the haline contraction coefficient (which ranges from 7.2 × 10 −4 to 7.9 × 10 −4 K −1 ) vary both spatially and temporally, depending on the ocean surface temperature and salinity. The remaining parameters are global constants: the acceleration due to gravity g, the reference salinity of seawater S ref , the reference specific heat of seawater c p , and the reference density of freshwater fw and seawater sw . Further discussion of the rationale behind our choice of the CORE forcing protocol and some discussion of relevant uncertainties associated with CORE appear in section 2.3.1. However, the formulation and sources of error arising from the bulk flux algorithms in the CESM model have been discussed extensively elsewhere (e.g., Yeager, 2004, 2009) and will not be addressed here.
There are only two small sources of error that are specific to the analysis of B in this manuscript, both of which arise from the fact that B is calculated offline. The first is due to neglected variations in coefficients and and associated correlations with variability in the heat and freshwater fluxes on submonthly time scales. Here, and are calculated offline from monthly mean sea surface temperature and salinity using the model equation of state (McDougall et al., 2003) and regridded in time to match the daily resolution of the heat and freshwater fluxes. The second source of error, which is only nonzero where sea ice forms and only impacts the subseasonal variability of B, is due to the absence of submonthly variability of the ocean internal heat/virtual salt gains due to frazil ice formation (see section 8.2 of Smith et al., 2010), which are only saved as monthly averages. We estimate that both of these sources of error are too small to impact the results.
In particular, we find that daily (vs. monthly) sea surface temperature variability and associated variability in , which is a far greater source of error than submonthly variability of , modifies the climatological mean buoyancy flux by less than 10% over 97% of the ocean surface and less than 1% over 75% of the ocean surface (not shown). In addition, the submonthly variations in modify the subseasonal variance of B by less than 10% over virtually the entire ocean surface and by less than 1% over 92% of the ocean surface (not shown). This error arising from submonthly variability of is generally greatest in major frontal zones, in which strong sea surface temperature gradients, vigorous mesoscale ocean variability, as well as large buoyancy fluxes of both signs all coexist.
Although we do not have the outputs of daily frazil ice formation to evaluate the effect of submonthly variability of frazil ice formation on subseasonal B variance, we find that frazil ice formation is mostly a very small contribution to the seasonal cycles of heat, freshwater, and buoyancy fluxes, even in regions that are seasonally covered by sea ice. For example, the maximum monthly mean heat flux associated with frazil formation is greater than 10 W/m 2 over only 2% of the ocean surface (not shown). In addition, the inclusion of frazil ice formation only changes the mean buoyancy flux by more than 10% over 3% of the ocean surface and by more than 1% over 9% of the ocean surface (not shown). So we do not expect submonthly variability in frazil ice formation to significantly impact the subseasonal variance of B either.

LP Experiment
The atmospheric state with smoothed winds is constructed in the frequency domain following the method for constructing the normal year (Large and Yeager, 2004). We apply a low-pass filter to pseudo momentum fluxes (also called pseudo-stresses) |u a |u a constructed from the zonal and meridional components of the normal year winds. The pseudo fluxes are then decomposed into harmonics with a fundamental frequency of 1 cycle per year via a Fourier transform. The amplitudes of all harmonics with a period shorter than 60 days are scaled down by a factor of 10 −2 . The resulting filtered pseudo fluxes are then inverse transformed, and the smoothed winds for the LP normal yearũ a are computed, for example, like the rescaled winds in the normal year (see section 9 of Large and Yeager, 2004). Throughout the paper, we use the term "smoothed winds" to refer toũ a , and "filtered winds" to refer to winds that have been filtered directly (without passing through pseudo fluxes) (e.g., as in Wu et al., 2016). The rationale behind our smoothing method is discussed in more detail in section 2.3.2.
The choice to use a filter cutoff time scale of 60 days to separate subseasonal from seasonal variability in the LP experiment is somewhat arbitrary. This 60-day time scale was chosen because it captures much of the subseasonal wind variability at all latitudes and leaves a reasonable number of harmonics (6) to represent the remaining long-period subannual wind variance. In the CORE normal year, the long-period (60-365 days) subannual harmonics represent the climatological annual cycle, which contains nonnegligible power in harmonics higher than the fundamental harmonic, as well as pure variability. In particular, it is important to note that some important modes of atmospheric variability, such as the Madden-Julian Oscillation (Zhang, 2005), project onto Harmonics 1-6 (i.e., 365-to 61-day periods), and there is variability at all time scales at all latitudes. Hence, the quantitative results are sensitive to the filter cutoff time scale.
It is impractical to quantify this sensitivity of the MLD to all possible cutoff time scales with the global mesoscale-resolving model. Instead, to provide some limited quantitative perspective on the sensitivity of the MLD to the filter cutoff time scale, results are presented below from several 1-D vertical column simulations (described in section 2.1.2) with a range of different cutoff time scales from 120 days to 1 day at a few locations. Future work is needed to quantify how each harmonic of subannual wind variability cumulatively impacts the MLD at all locations. However, many 1-D simulations, including some in areas with strong intraseasonal (30-90 days) atmospheric variability (not shown), indicate that the results are not highly sensitive to small changes in the cutoff time scale. Hence, we infer that the 3-D results are probably robust to small changes in the cutoff time scale and provide useful and general information about the impact of subseasonal winds on the MLD across the globe.

Rationale for the Experimental Protocol 2.3.1. CORE Forcing
The advantages and disadvantages of using an ocean-ice forcing data set like the CORE normal year for the CTL experiment have been discussed in detail elsewhere (Large and Yeager, 2004;Griffies et al., 2009). In addition, the CORE normal year forcing protocol has been tested in at least seven different ocean/sea ice models at a nominal 1 • resolution, and these solutions have been documented in the literature (Griffies et al., 2009) (see also Danabasoglu et al., 2014). Further, several high-resolution experiments with forcing similar to the CORE protocol on a nominal 0.1 • grid have been documented in close variants of the model configuration used here (Bryan and Bachman, 2015;Harrison et al., 2018;Maltrud et al., 2010) and in at least one other similarly high-resolution ocean model (Ibrayev et al., 2012).
This section briefly highlights some features of the CORE normal year that are particularly relevant for interpreting the results presented below. First, the forcing product does not contain interannual variability. Hence, the mean differences between the LP and CTL experiments can be quantified more precisely, because the differences are not obscured by the presence of interannual atmospheric variability. Second, the annual mean and seasonal cycle of the various components of the stress and buoyancy flux in the CORE normal year forcing products are reasonably well documented (Griffies et al., 2009;Large and Yeager, 2004), and the global patterns are, to a first approximation, consistent with other estimates, although there are significant outstanding uncertainties and differences, for example, in high-latitude buoyancy fluxes (Bourassa et al., 2013;Cerovečki et al., 2011). In addition, subseasonal wind variability has a significant rectified impact on the mean and seasonal cycle of the air-sea fluxes (e.g., Gulev, 1994;Hughes et al., 2012;Lin et al., 2018;Ponte and Rosen, 2004). Hence, it is important that this forcing product includes a climatological representation of temporal variability on periods ranging from 1 year to 12 hr (i.e., 6-hr resolution), and therefore, it includes the majority of the subseasonal atmospheric variance of interest to this study, including 2-to 10-day synoptic variability (wind speed power spectra are shown in the Figure S8).
The subseasonal variability of the fluxes of buoyancy and momentum are not as well understood compared to the annual mean and seasonal cycles, but they have been studied, both in the context of the NCEP reanalysis (e.g., Gulev and Belyaev, 2012;Zolina and Gulev, 2003) and in observations (e.g., Gille, 2005;Konda et al., 2010;Monahan, 2006;Monahan, 2008;Ogle et al., 2018;Sathiyamoorthy and Moore, 2002). For some analysis comparing the NCEP reanalysis winds that underpin CORE, in situ buoy winds, and satellite winds, see Figure S8 ( Atlas et al., 2011;Wentz et al., 2015), and see Figure 1 and section 3 of Alford (2001) (for further evaluation and discussion of the high-frequency winds at high latitudes, see also Alford, 2003). The upshot is that the subseasonal wind variance is fairly well represented in CORE, including the high-frequency wind variance at 0.5-to 5-day time scales. Additional validations at a few sites suggest that the CORE wind phases are coherent with observed satellite wind phases up to a frequency of about 1 cycle per day (not shown).
However, some small-scale, fast-moving mesoscale atmospheric variability is still missing from the forcing fields (the spatial resolution is ∼ 2 • ), which may result in reduced energy input from high-frequency winds to inertial oscillations (by perhaps as much as a factor of 3) and less near-inertial kinetic energy and mixing in the upper ocean with implications for the simulated ocean circulation and climate (Condron and Renfrew, 2013;Jochum et al., 2013;Li and Sriver, 2018;Rimac et al., 2013). In addition, since precipitation is only provided monthly, the subseasonal variability of the buoyancy flux is biased low in both experiments. Finally, the spatiotemporal relationships between the various atmospheric state variables are slightly corrupted in the normal year, since the amplitudes of the various Fourier harmonics for each state variable are changed to match climatological values over the Years 1958-2000, although the phase relationships are from the Year 1995 and are only modified near the start/end of the year to facilitate a smooth transition in multiyear simulations.

Filtered Pseudo Momentum Fluxes
In contrast to previous studies on the ocean's response to wind variability, in which the winds (e.g., Wu et al., 2016) or stresses (e.g., Lozier et al., 2008) were filtered, the winds are smoothed by filtering the atmospheric pseudo momentum fluxes |u a |u a (see section 2.2.2). There are two main advantages to forcing with winds that are smoothed by filtering the components of the pseudo momentum flux instead of forcing with filtered winds or filtered stresses. First, the annual mean pseudo momentum fluxes ⟨|u a |u a ⟩ t in LP are equal to those in CTL by construction. Hence, differences in the mean components of the stress at the air-ocean interface ⟨ ⟩ t [see (1)] between the CTL and LP solutions arise only from differences in the drag coefficient c d and the ocean currents u (⟨·⟩ denotes an average, and the subscript t indicates the averaging dimension is time). However, the dominant cause is presumably a reduction in the drag coefficient c d due to a reduction of the mean wind speed (e.g., Large and Pond, 1981;Large and Yeager, 2004), not a systematic increase in ocean current velocity in the direction of the wind ( Figure S5 shows that mean ocean surface kinetic energy is mostly reduced under LP). In addition, the online corrections to the wind magnitude and direction during the model runs do not introduce differences to ⟨|u a |u a ⟩ t between LP and CTL. Hence, the differences in ⟨ ⟩ t between LP and CTL are small. In particular, the components of ⟨ ⟩ t are reduced in magnitude by only about 10% or less in LP relative to CTL (Figures S1). Hence, mean properties that are linear in , such as the wind stress curl and Ekman suction, are modified by a similarly small percentage. The approximate equivalence of the time-mean components of the stress ⟨ ⟩ t in LP is beneficial because the large-scale ocean circulation depends strongly on ⟨ ⟩ t (more on this in the next paragraph), so it is informative to isolate the effect of subseasonal winds on ⟨ ⟩ t from the other effects of subseasonal winds on surface fluxes. Second, the sensitivity of the surface stress to ocean and ice velocity is retained, because the ocean/sea ice system is forced by atmospheric winds, and the fluxes are calculated online from bulk formulas. This is beneficial because this dependency significantly impacts the wind work on oceanic geostrophic flows, particularly the variable mesoscale flows, as well as other aspects of upper-ocean energetics and mixing (Duhaut and Straub, 2006;Dawe and Thompson, 2006;Eden and Dietze, 2009;Hughes and Wilson, 2008;Roberts et al., 2015;Renault, Molemaker, Gula, et al., 2016;. Several previous studies have demonstrated the strong impact of filtering the winds on the mean stress ⟨ ⟩ t , which is effectively unchanged in the LP experiment due to our different smoothing approach. For example, Munday and Zhai (2017) and Lin et al. (2018) show that, due to the nonlinear dependence of the stress (1) on the highly variable wind speed, about half of the zonal-mean zonal stress ⟨ x ⟩ tx is attributable to wind variability in the midlatitude storm tracks. Further, Wu et al. (2016) show that subseasonal winds drive substantial changes in the mean wind stress curl ∇×⟨ ⟩ t , which leads to significant changes in the barotropic ocean circulation as well as a roughly 50% increase in wind work on the geostrophic flow ⟨ ·u g ⟩ compared to a simulation where the submonthly winds are filtered directly (see also Zhai et al., 2012;Zhai, 2013). Munday and Zhai (2017) show that, as a result, the degree of wind variability modifies the sensitivity of the residual meridional overturning circulation to mean wind stress changes in an eddying channel that is an idealized representation of the Southern Ocean. In addition, Wu et al. (2016) find that this increase in the geostrophic wind work is primarily caused by increases in ⟨ ⟩ t · ⟨u g ⟩ t ; these increases are primarily caused by changes in the mean pseudo momentum fluxes, which are identical in the LP and CTL experiments presented here.
The motivating idea behind this study is that these CTL and LP experiments and future experiments forced by smoothed winds derived from filtered pseudo fluxes will serve to complement experiments forced by filtered winds and filtered fluxes. Together, these results will help to build a more complete understanding of the impacts of subseasonal winds on the ocean and thereby the role of subseasonal winds in coupled global Earth system dynamics. The present work, which documents the differences in the surface stress, surface buoyancy flux, and MLD between the LP and CTL experiments, is just a first step toward this broader goal.

MLD Definition and Validation Against Observations
There are many different definitions of the MLD. However, essentially all of the different definitions are a measure of the average vertical gradient of density or other hydrographic properties in the upper ocean. In addition, all of the definitions attempt to identify a sensible boundary between a relatively homogeneous near-surface layer and a deeper layer with a strong gradient (i.e., the pycnocline). Somewhat arbitrarily, we choose to define MLD using the definition of Large et al. (1997), which has been used in several recent studies of the mixed layer in high-resolution ocean models (e.g., Delman et al., 2018;Li and Lee, 2017;Jin et al., 2018). It is defined based on where the buoyancy frequency is maximum. First, the maximum buoyancy gradient relative to the surface (b k − b 1 )∕z k is obtained by searching over all discrete depths z k starting from the first grid cell which is at z 1 = −5 m. Then, the MLD is defined as the first depth where the local discrete buoyancy gradient, for example, (b k − b k−1 )∕(z k − z k−1 ), exceeds that maximum. We conduct most of our analysis using outputs of the daily maximum MLD to filter out diurnal variability, which is not well represented in these simulations because there is no diurnal cycle in the radiative forcing. The MLD from the 1-D simulations is calculated using the same definition as CESM, but the MLD is calculated offline from daily average temperature and salinity profiles.
Comparisons between the monthly climatology of the MLD in the CTL experiment and observations are presented in the supporting information for validation purposes ( Figure S3, see also de Boyer Montégut et al., 2004;Cabanes et al., 2013;McDougall and Barker, 2011;Peralta-Ferriz and Woodgate, 2015). Notably, the observational monthly MLD climatology is constructed using all individual Argo profiles from 2000 to 2017 so that both the observations and model use essentially the same MLD definition. The new Argo MLD climatology based on the Large et al. (1997) definition can be found online (see . Further, this new observational MLD climatology is compared with observational climatologies based on the same data but two other MLD definitions in Figure S7. The upshot is that the monthly climatology of the MLD is qualitatively well represented in the model and not especially sensitive to the MLD definition. However, there are still some quantitative biases relative to obervations, and the MLD does exhibit systematic quantitative sensitivity to the definition. Future work is needed to validate the simulated subseasonal variability of the MLD, but quantifying the subseasonal variability of the MLD with observations is beyond the scope of this study. In addition, future work is needed to quantify how the sensitivity of subseasonal MLD variability to subseasonal winds depends on the MLD definition. For example, we expect that MLDs defined with a very small density threshold (Brainerd and Gregg, 1995) or defined from vertical profiles of other tracers (e.g., bio-optical; see Carranza et al., 2018) may exhibit more variance at subseasonal time scales and/or greater sensitivity to subseasonal winds. However, we do not quantify the sensitivity of the model results to the choice of MLD definition because that would require rerunning the 3-D model.

Analysis Metrics
In section 3, several metrics are used to quantify the impact of subseasonal wind variability on the magnitude of the surface stress | |, the buoyancy flux B, and the MLD. The analysis code and model output are publicly available , and the method for computing these metrics is described briefly in this section.

Annual Mean
First, monthly means are computed from daily outputs at each model grid point. Then, in the case of B, the contribution due to frazil ice formation, which is only available as monthly output, is added. Then, the annual means for CTL are calculated at each grid point from the time-weighted arithmetic means of all 60 months, and the annual mean differences, CTL minus LP, are calculated from the time-weighted arithmetic means of all 60 monthly differences. Finally, a nonparametric bootstrap test with 1,000 samples is used to assess whether the annual mean differences are statistically different from zero (p < 0.05) at each grid point. Those grid points where the mean difference is not significantly different from zero are colored white in the plots. All the difference plots are presented as fractional differences, = ⟨CTL − LP⟩ t ∕⟨|CTL|⟩ t , that is the fractional contribution of subseasonal winds to the result in CTL (e.g., the fractional change in MLD is denoted MLD). The distributions of these standardized differences are not exactly normal, but the impact of using a nonparametric bootstrap approach rather than a t test is virtually indistinguishable in this analysis.

Seasonal Cycle 2.5.2.1. Amplitude
The seasonal cycle amplitude in each of the five simulated years in each simulation is calculated at each grid point from the 12 monthly means, for example, max t (B) − min t (B). Then the seasonal cycle amplitudes associated with each simulated year of the CTL simulation are averaged to produce a single estimate of the climatological seasonal cycle amplitude during the 5-year CTL simulation. Likewise, the five differences, CTL minus LP, are averaged to produce a single estimate of the climatological difference in the mean seasonal cycle amplitude between the CTL and LP simulations. As in the previous section, the bootstrap sampling approach is used to assess significance in the differences between the two simulations; this approach again yields results that are virtually indistinguishable from a t test.

Phase
The phase of the seasonal cycle is defined by the phase angle of the first harmonic of each 365-day annual time series (i.e., the harmonic with a 365-day period), which is derived from the Fourier transform of the annual time series during each of five years. The resulting phase angles from the five years are then processed like the seasonal cycle amplitudes to produce estimates of the climatological phase angle for the CTL and the climatological difference between the phase angle in CTL and LP at each model grid point.

Subannual and Subseasonal Variability
Unless otherwise specified, subannual variability indicates all variance with periods shorter than 365 days, subseasonal variability indicates all variance with periods shorter than 60 days, long-period subannual variability refers to variance with periods from 60 to 365 days that is not associated with the climatological seasonal cycle, and synoptic variability indicates all variance with 2-to 10-day periods. In addition, it is to be unstated but understood that since the output is daily averages, only periods from 2 days and up are explicitly included in the results. Yet high-frequency wind variability with periods from 0.5 to 2 days, which represents a small fraction of the subseasonal wind variance, is included in the forcing of the CTL experiment and removed from the forcing of the LP experiment. Hence, the indirect impacts of the high-frequency 0.5-to 2-day wind variability on longer time scales are included in the results to the extent that this high-frequency variability is fully represented in the CORE forcing (see sections 2.2.1 and 2.3.1). Finally, for all of the metrics described in this section, frazil ice formation is excluded from B (see section 2.2.1).
Both subannual and subseasonal variances (and standard deviations) are constructed in the frequency domain. First, five annual power spectra P( ) are created from the daily fields in each year. Then, five subannual variances are obtained for each grid cell in each experiment (LP and CTL) by integrating all frequencies in the power spectra except the zero-frequency (mean). The subseasonal variance includes only those frequencies higher than 1/60 cycles per day in the integral. The associated subannual and subseasonal standard deviations are computed by taking the square root of the respective integrals/variances. The results from the five years are then processed like the seasonal cycle amplitude to produce estimates of the climatological variances (and standard deviations) for the CTL and the climatological difference between the variances (and standard deviations) in CTL and LP at each model grid point. x + 2 calculated from the daily averaged components x and y during the last simulated year at subtropical (a,b) and subpolar (c,d) sites along 30 • W in the South Atlantic (SATL), as marked in Figures 13 and 18. The left panel shows results from the CTL (blue) and LP (red) 3-D simulations, whereas the right panel shows 30-day moving averages from eight 1-D simulations with different filter cutoff time scales, as indicated in the legend (in days). In addition, results are shown from the default CTL winds and the 3xHF winds, in which power in the pseudo-stress at all time scales shorter than 5 days is scaled up by a factor of 3.
The power spectra are also integrated in the frequency domain and normalized to obtain the cumulative fractional variance, that is where 1 =1/365 cycles per day is the first harmonic and n is the Nyquist frequency (1/2 cycles per day). This cumulative fractional variance CFV( ) quantifies the fraction of the subannual variance that is represented by frequencies lower than . For example, the fraction of subannual variance that is represented by subseasonal time scales is 1 − CFV( subs ) where subs = 1/60 cycles per day. And the fraction of subannual variance that is represented by synoptic time scales is 1 − CFV( syn ) where syn = 1/10 cycles per day. Hence, the fraction of subseasonal variance that is representated by synoptic time scales is (1 − CFV( syn ))∕(1 − CFV( subs )).

Results
The MLD is modified by ocean surface stress and buoyancy flux, which depend on the state of the atmosphere, ocean, and sea ice. Hence, we quantify the impacts of subseasonal winds on the annual mean, seasonal cycle, and subseasonal variability of the magnitude of the ocean surface stress | | (section 3.1), the ocean surface buoyancy flux B (section 3.2), and then the MLD (section 3.3). Each section begins with some examples to illustrate how subseasonal winds modify the annual cycles at two somewhat arbitrary sites in the subtropical and subpolar South Atlantic. Then the global results are reported.

Ocean Surface Stress 3.1.1. Local Examples
Time series comparing the magnitude of the surface stress | | in the CTL and LP experiments at two sites in the midlatitude South Atlantic highlight the significance of subseasonal wind variability for | | (Figure 1). The variance is larger, and the mean is greater at the subpolar site, but otherwise, the sites are qualitatively similar. At both sites, | | is dominated by subseasonal variability; hence, | | intermittently achieves both much larger and much smaller magnitudes in the CTL experiment compared to LP. In addition, due to nonlinearity, the 30-day moving average of | | is enhanced by up to about a factor of 2 in CTL compared to LP, even though the smoothing procedure holds the magnitude of the average stress |⟨ ⟩ t | (in contrast to ⟨| |⟩ t ) approximately constant. Further, 1-D simulations with various smoothing time scales applied to the wind show that all subseasonal time scales between 2 and 120 days generally enhance the 30-day moving average of | |. However, the synoptic time scales (shorter than 10 days) are particularly important at these  midlatitude sites: Scaling up the pseudo-stresses with periods shorter than 5 days by a factor of 3 has a positive effect on ⟨| |⟩ t that is greater in magnitude than the negative effect of smoothing all wind variability shorter than 120 days.

Global Annual Mean
Since the ocean surface stress is highly variable in time at most locations, the annual mean magnitude of the stress ⟨| |⟩ t is significantly reduced at most locations in the LP experiment compared to CTL (Figures 2a  and 2b). For example, the global average of ⟨| |⟩ t is reduced from 0.096 to 0.065 N/m 2 (32%) in LP relative to CTL. In addition, the percentage impacts are larger poleward of 30 • than at lower latitudes. In the middle-to-high latitude Northern Hemisphere (poleward of 30 • ), ⟨| |⟩ t is reduced by 50% (from 0.108 N/m 2 in CTL to 0.054 N/m 2 in LP), whereas in the middle-to-high latitude Southern Hemisphere, ⟨| |⟩ t is reduced by 40% (from 0.139 N/m 2 in CTL to 0.083 N/m 2 in LP). Further, in the Arctic and Antarctic (poleward of 66 • ), ⟨| |⟩ t is smaller in LP compared to CTL by 57% (0.043 N/m 2 ) and 52% (0.053 N/m 2 ), respectively.

Global Seasonal Cycle
In most locations, subseasonal (<60 days) variance of | | dominates the subannual (<365 days) variance in CTL (Figure 3). Hence, the seasonal cycle is not a dominant mode of temporal variability of | | in the global average power spectrum, and the first (365 days) harmonic explains only about 10% of the subannual variance on average (Figure 4f). Nevertheless, the amplitude of the seasonal cycle is comparable to or larger than the annual mean ⟨| |⟩ t in many locations (Figures 2a and 2c). In addition, the seasonal cycle and long-period subannual (60-365 days) variability is a dominant source of subannual variance of | | in a few locations, most notably the western tropical Indian Ocean and the far-eastern equatorial Pacific but also elsewhere equatorward of 15 • (Figure 3). Subseasonal winds have less of an impact on the seasonal cycle of | | compared to the annual mean ⟨| |⟩ t (Figures 2b and 2d). On average, the amplitude of the seasonal cycle is reduced from 0.10 N/m 2 in CTL to 0.09 N/m 2 in LP (10%). However, the response takes both signs (Figure 2d). The impacts are largest and most spatially coherent in the middle-to-high latitude Northern Hemisphere (poleward of 30 • ), where smoothing the winds reduces the amplitude of the seasonal cycle from 0.14 N/m 2 in CTL to 0.10 N/m 2 in LP (25%). Elsewhere, the impact of smoothed winds on the seasonal cycle in | | is relatively smaller and/or less spatially coherent.

Global Subseasonal Variability
Globally, about 75% of the subannual (<365 days) variance of | | is subseasonal (<60 days) in CTL, and this percentage is greater than 50% over 90% of the ocean surface area (Figure 3). The globally averaged power spectrum of | | and its cumulative integral show that synoptic time scales (<10 days) explain about 40% of the subannual variance, whereas the remaining subseasonal variance (10-to 60-day periods) explains about 35% of the subannual variance (Figures 4e and 4f).
Due to the dominance of subseasonal (<60 days) fraction of the subannual variance of | | in CTL, smoothing the subseasonal winds dramatically reduces the total subannual variance in LP compared to CTL (Figures 4b  and 4d). In fact, the global average subannual and subseasonal variances are reduced by 83% and 99%, respectively, in LP compared to CTL.

Ocean Surface Buoyancy Flux 3.2.1. Local Examples
At the two midlatitude South Atlantic sites, time series of the annual cycle of surface buoyancy flux B show that the seasonal cycle is the dominant mode of subannual variability, but subseasonal variability is also substantial ( Figure 5). In addition, the results show that subseasonal winds enhance the subseasonal variance of B. Yet the substantial subseasonal variability of B in the LP experiment also demonstrates that other subseasonal atmospheric and oceanic variability contributes significantly to the subseasonal variance of B, in contrast to the stress (cf. Figures 1 and 5). Various 1-D simulations with different smoothing time scales applied to the wind demonstrate that essentially all subseasonal wind variability produces a rectified increase in the amplitude of the seasonal cycle of B at these sites, but an impact on the annual mean ⟨B⟩ t is difficult to discern. At both the subtropical and subpolar sites, the qualitative effects of the subseasonal winds on B are similar, but the amplitudes of the seasonal cycle and variability of B are greater at the subtropical site, in contrast to the wind stress (cf. Figures 1 and 5).

Global Annual Mean
Subseasonal winds have a relatively small impact on the annual mean buoyancy fluxes ⟨B⟩ t compared to the annual mean stress ⟨| |⟩ t . The differences between LP and CTL are small (∼10%), take both signs, and are not significantly different from zero in many regions (Figure 6b). But some spatially coherent patterns emerge. For example, the buoyancy loss is reduced in LP relative to CTL in the vicinity of the midlatitude western boundary currents, where annual mean buoyancy loss is greatest in CTL. Buoyancy loss is also reduced in the Subantarctic zone of the Pacific sector of the Southern Ocean and deep convection sites in the subpolar North Atlantic. In contrast, much of the Arctic gains buoyancy on average as a result of subseasonal winds. In addition, subseasonal winds contribute to net buoyancy gain on the Grand Banks of Newfoundland and in the subpolar Atlantic and Indian sectors of the Southern Ocean.

Global Seasonal Cycle
The seasonal cycle is a dominant mode of temporal variability of the buoyancy flux over most of the ocean in CTL (Figure 7). Globally, the first subannual harmonic (a 365-day period) explains about half the total subannual variance of B on average (Figure 8f). In addition, the amplitude of the seasonal cycle of B is comparable to or larger than the annual mean ⟨B⟩ t , like | | (Figures 6a and 6c).
Subseasonal winds have a more modest impact on the amplitude of the seasonal cycle of B compared to | | (Figures 2d and 6d). Smoothing the winds reduces the global average amplitude of the seasonal cycle from 1.27 × 10 −7 m 2 /s 3 in CTL to 1.23 × 10 −7 m 2 /s 3 in LP (3%). However, the response takes both signs. At middle-to-high latitudes, subseasonal winds generally enhance the seasonal cycle amplitude by ∼10%, particularly in the subpolar North Atlantic and Pacific oceans and the far-western Pacific sector of the subpolar Southern Ocean (red in Figure 6d). For example, the seasonal cycle amplitude is reduced by 8% and 6% in LP relative to CTL poleward of 30 • in the Northern and Southern Hemispheres, respectively. Conversely, subseasonal winds suppress the seasonal cycle by ∼10% across most of the tropical Indian and Western Pacific oceans (±15 • latitude; blue in Figure 6d).

Global Subseasonal Variability
Globally, 38% of the subannual (<365 days) variance of B is subseasonal (<60 day) in CTL, and this percentage is greater than 50% over about 16% of the global ocean surface area (Figure 7). In addition, the subseasonal variability of B spatially mirrors the total subannual variability (Figures 8a and 8c). Further, B  Since subseasonal (<60 days) B variance is mostly less than longer-period (60-365 days) subannual variance, smoothing the winds on subseasonal time scales only reduces the globally averaged subannual B variance by 24% (Figures 8b and 8d). In addition, the globally integrated subseasonal (<60 day) variance of B is only reduced by 51% in LP relative to CTL (Figure 8d spatial variability (i.e., the shaded spread) of the MLD at these sites, except for a few months following the spring transition from deep to shallow MLDs.
The qualitatively similar increase in the MLD between the LP and CTL experiments in both the 1-D and 3-D simulations (dashed and solid lines, respectively, in Figures 9a and 9c) suggests that a substantial part of the impact of subseasonal winds on the MLD is due to changes in local surface buoyancy and momentum fluxes and 1-D vertical mixing physics. However, the differences (CTL-LP) are generally smaller between the 3-D experiments than between the 1-D experiments. Hence, 3-D processes substantially impact and generally suppress the sensitivity of the mean MLD to subseasonal winds at these sites.
Various 1-D simulations with different smoothing time scales demonstrate that all subannual time scales from 2 to 120 days have a discernible impact on the annual cycle of the MLD (Figures 9b and 9d). In particular, the mean MLD generally increases for each reduction in the filter time scale, from 120 to 2 days. In addition, the MLD increases substantially when the high-frequency (<5 days) pseudo-stresses are enhanced over the CTL in the 3xHF experiment. On the other hand, time scales longer than 30 days have a relatively small impact on the MLD compared to shorter time scales at these midlatitude South Atlantic sites, but the relative significance of these longer (30-120 days) time scales is greater at some other sites, such as the Labrador Sea deep mixing site and the western tropical Pacific where the Madden-Julian Oscillation is strong (not shown).

Global Annual Mean
Smoothing the subseasonal winds reduces the global annual mean MLD by 12.7 m from 66.7 m in CTL to 54.0 m in LP (19%) (red in Figure 10b). Although the mean MLD is almost everywhere deeper in CTL than LP, there is substantial spatial variability in the percentage difference. A significant fraction of the spatial variance in this difference between CTL and LP is explained by a zonally symmetric response, similar to the change in mean | | (Figures 2b and 10b). The change in the mean MLD is smallest in both absolute and percentage terms in the tropics and subtropics ( is that the global maximum is found in the Labrador Sea, where the annual mean MLDs are especially deep in CTL and therefore the absolute change in the MLD is maximal as well (Figures 10a and 10b).

Global Seasonal Cycle
The seasonal cycle is a dominant mode of temporal variability of the MLD across much of the ocean ( Figure 11). Globally, the first subannual harmonic (a 365-day period) explains about 60% the subannual variance of the MLD on average (Figure 12f). In addition, the amplitude of the seasonal cycle is greater than the annual mean on average (Figures 10a and 10c).
Smoothing the subseasonal winds reduces the amplitude of the seasonal cycle from 94.2 m in CTL to 80.5 m in LP (14%) on average (Figures 10c and 10d). In the tropics, subseasonal winds have a small impact on the seasonal cycle, and the western tropical Pacific and eastern tropical Indian ocean stand out as regions where the amplitude is actually larger by ∼10% in LP relative to CTL (similar to the buoyancy flux; see Figure 6d) (see also Waliser et al., 2003Waliser et al., , 2004. In the middle-to-high latitudes (poleward of 30 • ), smoothing the subseasonal winds reduces the amplitude of the seasonal cycle from 115.4 to 96.0 m (17%) in the Northern Hemisphere and from 138.9 to 109.8 m (21%) in the Southern Hemisphere (Figure 10d).
Example monthly climatologies of the MLD in different basins and at different latitudes reveal several qualitative differences between CTL and LP ( Figure 13): First, smoothing the wind generally reduces both the minimum (i.e., summer) and maximum (i.e., winter) MLD, but the absolute reductions are mostly greater when MLDs are deep (winter) and smaller when MLDs are shallow (summer), consistent with a reduced amplitude of the seasonal cycle, as in the South Atlantic examples. Second, the phase of the seasonal cycle is similar in LP and CTL. More generally, a map of the phase shift between CTL and LP associated with the first harmonic (with a 365-day period) shows that subseasonal winds generally have a modest impact on the phase of the seasonal cycle ( Figure 14). Notable regions where subseasonal winds delay the time when the MLD reaches its deepest depth include the subpolar North Pacific and the Antarctic marginal ice zone.
On a pointwise basis, the changes in the climatological seasonal cycle between LP and CTL are mostly statistically insignificant in the regional hot spots of deep winter MLDs (which are defined here by the top 7% of grid points in terms of annual mean MLD; Figure 10d). However, suppressing the effects of intrinsic oceanic variability by averaging all these grid points into regional hot spots shows that the seasonal cycle of MLD is reduced in the LP scenario (relative to CTL) in all these regions ( Figure 15). The Labrador Sea is most notable (as in Holdsworth and Myers, 2015). There, the seasonal cycle amplitude exceeds 1,000 m in some locations, and the subseasonal winds are crucial to sustaining this large seasonal cycle amplitude in CTL. At other hot spots of deep winter MLDs besides the Labrador Sea, the effect of subseasonal winds is smaller, but at each hot spot, subseasonal winds contribute to deeper winter MLDs and a larger seasonal cycle amplitude on average. To a first approximation, the phase of the climatological seasonal cycle is not strongly modified by subseasonal winds at these hot spots with deep winter MLDs. However, the Northern Hemisphere subtropical mode water formation site near the Kuroshio is notable in that subseasonal winds cause the mixed layer to shoal slightly later in the spring (in CTL compared to LP).

Subseasonal Variability
Globally, subseasonal (<60 days) MLD variance represents only 17% of the subannual (<365 day) MLD variance in CTL, and subseasonal MLD variance represents more than 50% of the subannual variance of the MLD over only 5% of the global ocean surface area (see Figure 11 and contrast with Figures 3 and 7). The regions of relatively high subseasonal MLD variance are loosely concentrated in regions of high eddy kinetic energy and/or deep winter MLDs and relatively high subseasonal B variance (cf. Figures 11 and 12c with Figures S5, 10c, and 7). The equatorial oceans are a prime example.
Smoothing the subseasonal wind only reduces the globally averaged subseasonal MLD variance by 26% and the subannual variance by 27% (not shown; changes in the associated standard deviations are plotted in Figure 12d). Both the fractional and absolute impacts of smoothing the winds are greater at middle-to-high latitudes, but significant subseasonal variability in MLD remains at all locations. In other words, about three fourths of the global subseasonal MLD variability is explained by factors other than subseasonal wind variability, such as intrinsic subseasonal ocean/sea ice variability and subseasonal variability of the atmospheric surface temperature and humidity. The globally averaged power spectrum and its cumulative integral (Figures 12e and 12f) suggest that the impacts of subseasonal wind variability on global MLD variability is relatively small at all subseasonal time scales, but the impacts are slightly stronger at synoptic time scales.
On a pointwise basis, the changes in subseasonal variability in MLD are not statistically significant in the regional hot spots, except for the Labrador Sea. However, by calculating the change in the MLD power spectra averaged over the regional hot spots, more robust results can be obtained (Figure 16). The Labrador Sea stands out starkly in these comparisons, as a region where MLD variance is damped by about a factor of 5 Figure 13. Monthly climatologies of the MLD averaged over several 3 • by 5 • boxes located in the (a,c,e,g,j,m,p) Atlantic, (h,k,n,q) Indian, and (b,d,f,i,l,o,r) Pacific ocean basins as marked on the map. The control run (CTL) is marked by black line, the low-pass run (LP) by dashed, and the difference (CTL-LP) in red, for which ±2 standard errors of the mean are shown (assuming five independent estimates of the mean in each box, one from each year). Where the red lines do not bracket zero, the difference CTL-LP is statistically significant. For comparison, gray lines are from an Argo climatology, which is derived using the same MLD definition as the model (see supporting information). across all periods from 2 to 365 days. In addition, there is a fat tail at low spectral power in the Labrador Sea of the LP simulation: Some grid cells have especially low variance without subseasonal winds. We hypothesize that the fat tails arise because deep winter mixing exhibits a threshold-like (on/off) behavior (not shown), and the presence or absence of deep winter MLDs has a dramatic impact on subseasonal MLD variability. But future work is needed to fully investigate how the winter MLDs respond to subseasonal winds in the Labrador Sea, and a thorough analysis of this region is beyond the scope of this paper. In contrast, all the other hot spots exhibit very slightly or negligibly smaller subseasonal variance in LP compared to CTL.

Discussion
This section discusses the possible causes of the sensitivity of both the climatological mean and subseasonal variability of the MLD to subseasonal winds. In addition, areas of uncertainty that are in need of attention from more narrowly focused process studies are highlighted.

Mean MLD Response
Subseasonal winds deepen the mean MLD, as shown in Figure 10b, in several ways. The following sections consider possible causes of this rectified effect.

Changes in the Surface Stress
As shown in section 3.1, subseasonal winds contribute to the average magnitude of the ocean surface stress ⟨| |⟩ t as well as the subseasonal variability of | |. First, we discuss the changes in the mean and then the variability of | |.

Figure 15.
Monthly climatologies of the MLD averaged in regional hot spots (a) of deep mixing (top 7% in annual mean MLD), the same regions shown in Figure 16. These regions are found in the subantarctic zone of the Southern Ocean (b), south of the Kuroshio (c), the Labrador Sea (d), south of the Gulf Stream (e), and the Greenland-Iceland-Norwegian Seas (f). The control run (CTL) is marked by black line, the low-pass run (LP) by dashed, and the difference (CTL-LP) in red, for which ±2 standard errors of the mean are shown (assuming five independent estimates of the mean in each hot spot, one from each year). Where the red lines do not bracket zero, the difference CTL-LP is statistically significant. For comparison, gray lines are from an Argo climatology, which is derived using the same MLD definition as the model (see supporting information).

Figure 16.
Power spectra of MLD averaged in regional hot spots of deep MLDs, defined as in Figure 15. Area-weighted means are indicated by solid lines in (b)-(f); the spread between the 5th and 95th percentiles are shaded.
By enhancing the mean magnitude of the stress ⟨| |⟩ t , subseasonal winds increase the average mechanical energy available to mixed layer turbulence by strengthening the mean shear in the upper ocean, including under sea ice (for a review, see Large et al. (1994)), and thereby, subseasonal winds can deepen the MLD. In particular, theoretical dimensional and scaling arguments as well as idealized process simulations and laboratory experiments suggest that the MLD may scale with the mean friction velocity u * = √ | |∕ sw where sw is the density of seawater (e.g., Kato and Phillips, 1969;Price, 1998;Price and Sundermeyer, 1999;Zilitinkevich et al., 2002;McWilliams et al., 2009;Mironov and Fedorovich, 2010;Yoshikawa, 2015;Pham and Sarkar, 2017). For example, scaling and dimensional arguments alone suggest that the MLD scales with u * ∕f in an idealized unstratified rotating fluid column, the MLD scales with u * ∕N in a nonrotating-stratified fluid column, and the MLD scales with some combination of u * ∕f and u * ∕N (e.g., u * ∕ √ N; (Pollard et al., 1972;McWilliams et al., 2009)) in a rotating-stratified fluid under a neutral buoyancy flux B = 0 (f is the Coriolis frequency, and N is the buoyancy frequency). In addition, under very strong stable buoyancy forcing B ≫ 0 and relatively small f and N, the sensitivity of the MLD to u * may be even greater, since the Monin-Obukhov depth scale u 3 * ∕B (Monin and Obukhov, 1954) may be the relevant dimensional depth scale in this limit. In contrast, in regions of deep MLD under strong surface buoyancy loss B ≪ 0 (such that |u 3 * ∕B| << MLD), scaling arguments suggest that u * has little influence on the MLD. In order to identify the impacts of subseasonal winds on MLD that are attributable to changes in the mean u * , we compare the fractional changes u * and MLD, where = ⟨CTL−LP⟩ t ∕⟨|CTL|⟩ t , in the 3-D simulations. We find that the reduction in the global area and annual averaged friction velocity of 15% from 0.0088 m/s in CTL to 0.0075 m/s in LP is comparable to the corresponding 19% reduction in the MLD (cf. Figures 10  and S6). In addition, about 20% of the spatial variance in the response of the climatological annual mean MLD to subseasonal winds is explained by the change in mean u * (Figure 17).
However, it is notable that the remaining variance in MLD− u * has some coherent large-scale spatial structure ( Figure 18). In both hemispheres, the subpolar and equatorial latitudes (45-60 • and <15 • , red in Figure 18) are regions where the MLD is more sensitive to subseasonal winds than u * . Conversely, subtropical and polar latitudes (∼30 • and >60 • , blue in Figure 18) are regions where u * is generally more sensitive to subseasonal wind than the MLD, although the Antarctic seasonal ice zone is not entirely consistent with this zonally averaged perspective.
The near zonal symmetry of MLD − u * suggests a relatively simple explanation for these differences. For example, one possible cause of the zonal bands of lower MLD − u * near 30 • is the relatively large seasonal and subseasonal variability of the surface buoyancy flux, which we discuss further in the next section. In addition, subseasonal winds are associated with a positive mean buoyancy flux in the Arctic, which may lower MLD − u * there. However, the departures from the relationship u * ∼ MLD are not readily attributed to changes in stratification between LP and CTL; the fractional differences in u * ∕N, u * ∕f, and u * ∕ √ N produce results similar to those of u * (not shown).
Subseasonal variability of the stress or friction velocity, which is almost entirely caused by subseasonal winds (Figure 4b), may also have a rectified effect on the mean MLD. McWilliams et al. (2009) find that the rectified effect of longer-period synoptic variability in the wind speed (up to 10-day time scales) is a modestly shallower mean MLD. Hence, we hypothesize that 5-to 60-day stress variability acts to suppress the sensitivity MLD. On the other hand, high-frequency winds with frequencies near the Coriolis frequency f are resonant with ocean/sea ice inertial oscillations and therefore inject a relatively large amount of kinetic energy into the upper-ocean/sea ice system (compared to their fraction of the wind variance) (Alford, 2003;Alford et al., 2016;Jochum et al., 2013;Pollard and Millard, 1970;Rimac et al., 2013;Roberts et al., 2015), and  Figure 17), that is the magnitude and sign of the vertical distance from the 1-to-1 line in Figure 17b at every location. At red points, the fractional change MLD is greater than the fractional change u * , and the point is above the 1-to-1 line in Figure 17b, but at blue points, the fractional change MLD is less than the fractional change u * , and the point is below the 1-to-1 line in Figure 17b. The physical locations of the two 1-D vertical column simulations are marked as magenta circles for reference. these high-frequency winds deepen the mean MLD in idealized 1-D column simulations . However, the results from the 1-D simulations presented above suggest that removing high-frequency wind variability on periods shorter than 5 days has relatively little impact on the mean MLD (Figures 9b and  9d). Yet the 1-D experiments also show that amplifying the high-frequency wind variability can increase the mean MLD, as in McWilliams et al. (2009). Hence, future work is needed to more fully assess the global sensitivity of the mean MLD to high-frequency winds, but high-frequency winds in the CORE forcing contribute relatively little to MLD in our global simulations.

Changes in the Buoyancy Flux
The sensitivity of the mean buoyancy flux to subseasonal winds B is both weak and not obviously related to MLD (cf. Figures 6b and 10b). Some notable exceptions, where B may have a significant impact on MLD, include parts of the Arctic, where B is positive and MLD− u * is relatively low, and the Labrador Sea, where B is negative and MLD− u * is relatively large.
On the other hand, the amplitude of the seasonal cycle of B is enhanced by subseasonal winds, particularly at middle-to-high latitudes where MLD is also largest (Figures 6d). In addition, subseasonal B variance is enhanced by subseasonal winds at virtually all locations, similar to MLD, but the fractional increase in subseasonal B variance is greater at low latitudes than at midlatitudes (Figure 8d), in contrast to MLD. In addition, theoretical (Garwood, 1977;Kraus and Turner, 1967) and idealized 1-D numerical studies  suggest that increases in the seasonal and subseasonal B variances are associated with a shallower mean MLD. More generally, seasonal and subseasonal B variances (e.g., in the subtropics) may reduce the sensitivity MLD− u * . This would imply that the increase MLD due to subseasonal winds in the 3-D simulations is smaller than what would be obtained if the mean u * were increased in idealized process studies without B variability. However, future work is needed to clarify the impact of B variability on MLD− u * .

3-D Ocean Dynamics
Another factor that may influence the sensitivity of the mean MLD to subseasonal winds is 3-D ocean dynamics that is not captured by the scaling arguments and 1-D ideas discussed above but resolved in the 3-D CTL and LP simulations. For example, by modifying the energetics of mesoscale eddies (e.g., Barton et al., 1993;Duhaut and Straub, 2006) (see Figure S5), subseasonal winds can modify the mixed layer restratification process (Gelderloos et al., 2011;Katsman et al., 2004) and thereby the mean MLD. In addition, subseasonal wind variability may have a rectified impact on the mean ocean circulation and thereby the mean stratification and MLD, even though the mean zonal and meridional stress ⟨ ⟩ t and therefore the large-scale Ekman pumping/suction are approximately the same in the LP and CTL simulations.
Differences between 1-D and 3-D simulations of the same location represent one measure of how 3-D ocean dynamics influences the sensitivity of the MLD to subseasonal winds. For example, the two midlatitude sites in the South Atlantic (shown in Figures 1, 5, 9, and 19) represent two different regimes in the relationship between u * and MLD (Figure 18). In the subtropical South Atlantic, MLD is very close to u * in the 3-D experiments, and both the mean and seasonal variability of MLD are fairly well captured by u * (Figure 19a). However, the magnitude of MLD is more than 50% greater than u * in the analog 1-D experiments. In the subpolar South Atlantic, both the 3-D and 1-D models suggest that MLD is greater than u * , but again, the sensitivity MLD is much greater in the 1-D model than in the 3-D model (Figure 19b). These results demonstrate that, despite theoretical scaling arguments that u * ∼ MLD, the global correlation in a mesoscale-resolving ocean model is nontrivial to interpret. 3-D ocean dynamics evidently have an important influence on the sensitivity of the climatological mean MLD to subseasonal winds. In particular, these results suggest that 3-D ocean dynamics reduce the sensitivity MLD− u * (i.e., they reduce the slope of the regression lines in Figure 17). In addition, 3-D ocean dynamics may even enhance the correlation u * ∼ MLD (i.e., reduce the range of different MLD values associated with a given u * , as in Figure 19a).
Future work is needed to more thoroughly understand how 3-D ocean dynamics influences this sensitivity of the MLD to subseasonal winds and to generalize these results beyond a few sites. But it is notable that these results are reminiscent of observations that the MLD and u * ∕ √ N are correlated on coastal continental shelves (Lentz, 1992), where surface buoyancy fluxes and 3-D advection may be important as well. Hence, some existing observational evidence supports the hypothesis emerging from the model that MLD ∼ u * , even when advection and buoyancy fluxes are also important.

Subseasonal MLD Response
Subseasonal variability represents only 17% of the total subannual variance of the MLD, and subseasonal winds are responsible for less than a quarter of this subseasonal MLD variance. This small contribution of the subseasonal time scales to the total subannual variance of MLD (17%) compared to the surface forcing by | | (75%) and B (38%) presumably reflects the fact that the MLD is set by the integrated effects of | | and B. That is, the surface-forced subseasonal MLD variability in CTL is expected to have a frequency power spectrum proportional to −2 or −3 , given that the subseasonal B and | | power spectra are proportional to 0 to −1 and therefore a lower fraction of the subannual variance at subseasonal time scales (e.g., Hasselmann, 1976). However, although subseasonal surface stress variance is effectively eliminated in LP and subseasonal buoyancy flux variance is reduced by about 50%, about three fourths of the subseasonal variance of the MLD in CTL remains in LP (Figures 4c, 4d, 12c, and 12d). These results suggest that there is a nonlinear relationship between MLD and the surface forcing that is not captured by the linear stochastic model of Hasselmann (1976) and/or subseasonal variance of the MLD is more strongly controlled by factors other than subseasonal surface fluxes of momentum or buoyancy, such as intrinsic mesoscale oceanic variability. In addition, it may be noted that the smallness of subseasonal MLD variability compared to the mean and seasonal cycle supports an important role for ocean dynamics, since the linear stochastic model dynamics are expected to be a better representation of subseasonal MLD variability where the subseasonal variability is a small perturbation to the climatological seasonal cycle.
Future work is needed to explore the hypothesis that 3-D ocean dynamics is the dominant cause of global subseasonal MLD variability. But the comparisons between the 3-D and 1-D simulations in the South Atlantic (Figure 9 and section 3.3.1) support the hypothesis that subseasonal MLD variability is strongly controlled by ocean dynamics in these sites. In particular, the 3-D simulations show that the mesoscale spatial variability (represented by the shading in Figures 9a and 9c) is small compared to the amplitude of the seasonal cycle but comparable in magnitude to the subseasonal temporal variability at a single grid point (subseasonal wiggles in the solid lines). In addition, neither the mesoscale spatial variability nor the subseasonal variability at a single grid point are discernibly reduced (relative to CTL) when removing subseasonal winds in LP, except perhaps near the spring transition. Finally, the subseasonal variability of the MLD is also discernibly smaller in both of the 1-D simulations (CTL or LP) compared to either analog 3-D simulation, except perhaps during the spring. At both South Atlantic sites, these results are consistent with the hypothesis that 3-D ocean dynamics is the dominant driver of subseasonal MLD variability and subseasonal wind variability has comparatively little, albeit nonnegligible, impact on subseasonal MLD variability.
Future work should also explore how the MLD response to subseasonal winds depends on spatial scale as well as time scale. Our hypothesis is that the subseasonal variability of the MLD may be driven by different dynamics at different spatial and temporal scales that have not been separated by the above analysis, which focuses on the finest resolved spatial scales (3-11 km). One idea is that the ocean controls subseasonal MLD variability more strongly on scales corresponding with mesoscale eddies (100 km, 10-60 days) where mesoscale ocean turbulence is energetic, and the atmosphere controls subseasonal MLD variability more strongly on scales corresponding with synoptic storms (1,000 km, 2-10 days) where synoptic atmospheric activity is vigorous. Another idea is that subseasonal anomalies in surface stress may control MLD variability more strongly on shorter synoptic time scales (2-10 days), but subseasonal anomalies in buoyancy flux may be a more important driver of subseasonal variability in MLD on longer time scales (10-60 days). The idea of scale-dependent control is partially motivated by an analysis showing scale-dependent oceanic control of sea surface temperature and surface heat flux by Bishop et al. (2017), although MLD is not trivially related to sea surface temperature and/or heat flux.

Poorly Represented and Excluded Processes
Processes that are poorly represented or excluded in these simulations may cause the sensitivity of the MLD to subseasonal winds to deviate from the results reported above in the real ocean or other simulations. For example, subseasonal winds may modify Langmuir (Fan and Griffies, 2014;Li et al., 2016;Van Roekel et al., 2012) and submesoscale (Fox-Kemper et al., 2011;Thomas et al., 2013Thomas et al., , 2016, 2019 turbulence, neither of which are represented in the models used in this study. In addition, subseasonal winds modify inertial oscillations and internal waves, which are only marginally represented in the global model (see sections 2.2.2 and 2.3.2). Finally, as discussed in section 2, some features of the real atmosphere, such as intense mesoscale cyclones and fronts (<200 km) and submonthly precipitation variability are not included in CORE, and these missing features of the subseasonal forcing may systematically modify the MLD and its sensitivity to subseasonal winds. Hence, future work is also needed to determine the significance of the poorly represented or excluded effects for the impacts of subseasonal winds on the MLD.

Summary and Conclusions
It is well known that subseasonal variability of the surface wind strongly influences atmospheric surface fluxes of heat, freshwater, and momentum (e.g., Wanninkhof, 1992;Gulev, 1994;Zolina and Gulev, 2003;Ponte and Rosen, 2004); however, the impacts of subseasonal winds on the ocean are not fully understood. Motivated by the hypothesis that subseasonal surface wind variability has a significant but poorly understood impact on ocean surface MLDs with implications for climate dynamics and biogeochemistry, this paper quantifies the sensitivity of the climatological mean, seasonal cycle, and subseasonal variability of the ocean surface stress, buoyancy flux, and MLD to subseasonal surface wind variability in a global mesoscale-resolving ocean/sea ice numerical circulation model for the first time. Complimentary results from analog 1-D vertical column simulations at two specific sites in the midlatitude South Atlantic help to distinguish 1-D and 3-D controls on the MLD as well as the different time scales influencing the sensitivity of the MLD to subseasonal winds.
To our knowledge, no previous modeling study has separated the effects of subseasonal winds from the ocean surface forcing by filtering the pseudo-stresses, as we have done here. The novelty of our approach is that it allows us to approximately isolate the effects of subseasonal winds on the magnitude of the stress ⟨| |⟩ t from their effects on the mean zonal and meridional components of the stress ⟨ ⟩ t and retain surface flux feedbacks. Hence, confounding changes in ocean circulation due to changes in ⟨ ⟩ t are reduced compared to experiments with ocean models forced by filtered winds (e.g., Wu et al., 2016), and some air-ice-ocean interactions are retained that are missing from ocean models forced by filtered surface fluxes (e.g., Lozier et al., 2008). However, the main conclusions of this paper are not expected to be especially sensitive to our choice of smoothing procedure (comparisons with Wu et al. (2016) support this idea).
The main conclusion of the paper is that the mean MLD is relatively sensitive to subseasonal winds to the extent that they modify the mean friction velocity u * = √ | |∕ sw , but the subseasonal MLD variance is relatively insensitive to subseasonal winds despite their strong impact on local subseasonal variability of the surface momentum and buoyancy fluxes. In particular, we find that smoothing the subseasonal wind variability shorter than 60 days reduces the globally averaged annual mean MLD by about 19% or 12.7 m. This reduction in the mean MLD is closely related to the corresponding (15%) reduction in the mean u * . In fact, about a fifth of the spatial variance in the response of the mean MLD to subseasonal wind, including the relatively larger response at middle-to-high latitudes compared to the tropics, can be explained by the response of the mean u * . Future work is needed to more fully understand and observe how 3-D ocean dynamics as well as subseasonal variability of the buoyancy and momentum fluxes modify the effect of a stronger mean u * to produce a deeper mean MLD with subseasonal winds. But the results presented above and existing theories support the hypothesis that all of these other factors tend to counteract or suppress the effect of the enhanced mean u * and reduce the sensitivity of the mean MLD to subseasonal winds.
We also find that subseasonal MLD variance is weak; it only represents about 17% of the total subannual MLD variance, which is dominated by the seasonal cycle ( Figure 11). In addition, although the subseasonal winds cause virtually all of the subseasonal variability of the surface stress and half of the subseasonal variability of the surface buoyancy flux, we find that subseasonal winds are only responsible for about a quarter of the subseasonal MLD variance (Figure 12). Hence, we conclude that the subseasonal MLD variance is set primarily by intrinsic ocean dynamical variability at the spatial scales considered here (3-11 km), particularly in regions of high mesoscale kinetic energy.

10.1029/2019JC015166
In addition, we find that both the amplitude and phase of the seasonal cycle of the MLD can be significantly modified by subseasonal winds. Subseasonal winds mostly increase the amplitude of the seasonal cycle of the MLD at middle-to-high latitudes, but the western tropical Pacific and eastern tropical Indian oceans are notable regions where subseasonal winds suppress it. The modifications to the phase of the seasonal cycle of MLD is mostly less than 1 month or statistically insignificant, but there are a few significant exceptions such as the subpolar North Pacific where the deepest MLDs occur later in the year due to subseasonal winds.
Finally, previous modeling studies have shown that the deep winter MLDs in the Labrador Sea are highly sensitive to subseasonal wind forcing (Holdsworth and Myers, 2015;Wu et al., 2016). We also find that deep mixing in the Labrador Sea is highly sensitive to subseasonal winds, but we also show that other sites of deep MLDs around the globe are not especially sensitive to subseasonal winds. That is, the Labrador Sea is relatively unique in its sensitivity to subseasonal wind variability and deserves further study.
More comprehensive analysis of simulations and observations is needed to fully understand how atmosphere, sea ice, and ocean dynamics interact to set the mean and variability of the MLD. However, these results provide global and process-aggregated perspective to guide future studies. For example, the results suggest that detecting direct relationships between subseasonal wind and MLD variability is difficult, particularly at a single location (e.g., from a mooring). Hence, future studies of these direct effects should be targeted at specific regions where the impacts of subseasonal winds on subseasonal MLD variability are most readily detectable in the simulation results above (e.g., in regions where subseasonal MLD variance is strong and highly sensitivity to subseasonal winds, but ocean eddy kinetic energy is weak) and plan for a small signal to noise ratio. In addition, the results suggest that future work should more fully explore the interannual relationships between subseasonal wind variance and the seasonally averaged friction velocity and MLD in both observations and models. Finally, the results demonstrate that the relationships between subseasonal atmospheric variability and the mean MLD are sensitive to 3-D and possibly even mesoscale ocean dynamics. Hence, future work should evaluate the degree to which mesoscale ocean dynamics modulates the sensitivity of the mean MLD to changes in atmospheric forcing and ensure that these sensitivities are accurately represented in models that do not resolve ocean mesoscales. the Climate and Global Dynamics Laboratory, made this collaboration possible. The CESM project is supported primarily by the NSF. This material is based upon work supported by NCAR, which is a major facility sponsored by the NSF under Cooperative Agreement 1852977. Computing and data storage resources, including the Cheyenne supercomputer (doi: 10.5065/ D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. We thank all the scientists, software engineers, and administrators who contributed to the development of CESM2. Finally, this work would not have been possible without the hard work of many who develop the CESM. CCMP Version-2.0 vector wind analyses are produced by Remote Sensing Systems. Data are available online (www.remss.com). Argo data are available at http://marine.copernicus.eu/ services-portfolio/access-to-products/ (Cabanes et al., 2013).