Submesoscales Enhance Storm‐Driven Vertical Mixing ofNutrients: Insights From a Biogeochemical Large EddySimulation

Introduction

In the midlatitudes, strong winds associated with storms intermittently enhance upper-ocean turbulence and hence vertical mixing of nutrients and phytoplankton (e.g., Carranza & Gille, 2015;Marra et al., 1990;Nicholson et al., 2016;Rumyantseva et al., 2015). These storm events can cause significant intraseasonal variability of upper-ocean phytoplankton concentrations and therefore significant divergence from the mean seasonal cycle (e.g., Thomalla et al., 2011;Waniek, 2003). These variations in phytoplankton phenology impact the higher trophic levels of marine ecosystems (e.g., Cushing, 1990;Platt et al., 2003).

An unresolved question is to what degree quasi two-dimensional ocean submesoscales (defined here descriptively to refer to features with horizontal length scales from 0.1 to 10 km, and small aspect ratio) interact with smaller-scale wind-driven three-dimensional upper-ocean turbulence (<100 m) to modify the physical-biogeochemical response to a storm. This study addresses this question in the oceanic context of a fairly typical autumn storm that was observed in the midlatitude North Atlantic at 48.78 • N, 16.38 • W during the OSMOSIS project (e.g., Lucas et al., 2019;Painter et al., 2016;Rumyantseva et al., 2015) by showing how submesoscales impact net (submesoscale + turbulent) vertical nutrient flux and resulting biological production during and shortly after the passage of the storm in an ocean model.

In the absence of storms, prior work has shown that submesoscales can both enhance vertical tracer fluxes (e.g., Capet et al., 2008;Lévy et al., 2001;Mahadevan & Archer, 2000) or suppress vertical fluxes by restratifying the ocean mixed layer or reducing the depth of the surface mixing layer in ocean models (e.g., Bachman et al., 2017;Boccaletti et al., 2007;Fox-Kemper et al., 2008;Mahadevan et al., 2012). Enhanced vertical fluxes, for example, can induce a net increase in the nutrient delivery to the well-lit euphotic layer These conclusions are essentially based on regional and process-study primitive equation ocean models, in which submesoscales are at least partially resolved but small-scale turbulence is unresolved and parameterized. It is not clear whether submesoscale interactions with turbulent scales, and therefore the net impacts of submesoscales on biogeochemistry during storms, are realistically represented in such models.

Large eddy simulations (LESs) that explicitly resolve some submesoscales as well as smaller-scale turbulence have started providing additional insights into the interactions between submesoscales and turbulence. For example, such studies have shown energy transfers from submesoscales to small turbulent scales via submesoscale instabilities (e.g., Bachman et al., 2017;Skyllingstad & Samelson, 2012;Skyllingstad et al., 2017;Thomas & Taylor, 2010). If turbulent scale motions are more energetic as a result, this may enhance vertical fluxes. In a recent paper, Whitt and Taylor (2017) used LES to explore the physical response of submesoscales coupled with small-scale turbulence to a storm. Their simulations showed that energetic wind-driven submesoscale motions nonuniformly restratified the mixed layer while deepening the reach of boundary-layer turbulence during the storm. This result ensued from spatial heterogeneity in stratification and turbulence.

Here, we build on this previous study (Whitt & Taylor, 2017) by describing the influence of submesoscales coupled with small-scale turbulence on biogeochemistry during and shortly after a strong storm using biogeochemical LES (BLES; for some other recent examples of BLES without submesoscales, see Lewis et al., 2017or Brereton et al., 2018; for an LES with submesoscales and a plankton model but without nutrients or wind, see Taylor, 2016; or for submesoscales and passive tracers under weak atmospheric forcing, see Smith et al., 2016). Results are obtained by comparing simulations in which submesoscales are present or absent. These results highlight how submesocale and smaller-scale turbulent nutrient fluxes are modified by the interaction between submesoscale and smaller-scale turbulence, and how this impacts biological production and nutrient and phytoplankton vertical and horizontal distributions. In a second step, we evaluate the efficacy of existing commonly used parameterizations of submesoscales (Fox-Kemper et al., 2008) and small-scale turbulence (Large et al., 1994) in a column model against the BLES benchmarks. We find that the nutrient flux by small-scale turbulence is enhanced by submesoscale interactions during the storm but suppressed shortly after the storm, with the consequence of an increased net biological production, and that the commonly used parameterizations fail at capturing these effects. These results motivate future physical/biogeochemical theory and parameterization development, which can be validated by comparison with these BLES benchmarks, and provide a framework for interpreting the simulation results, as well as additional BLES simulations to more fully quantify the sensitivity of the simulation results to the parameters that define the numerical experiments.

Models 2.1.1. BLES Setup

The BLES resolves the majority of the variance associated with three-dimensional fluid motions and is used to evaluate the impact of submesoscale/turbulent-scale interactions on biogeochemestry. The rotating stratified Boussinesq equations are solved on an f -plane without surface wave effects over horizontal wavelengths ranging from about 2 km to about 4 m, which include both submesoscales and smaller-scale turbulence. The numerical model, spatial domain, and physical initial condition (of density, velocity, and pressure) are as described in Whitt and Taylor (2017; see also Taylor, 2008). Briefly, the model solves the Boussinesq equations together with the evolution equations for the four reactive biogeochemical tracers using a low-storage third-order Runge-Kutta time marching method, a pseudospectral method to represent horizontal derivatives, and second-order centered differences to represent vertical derivatives. Fluxes of momentum associated with subgrid-scale turbulence with wavelengths smaller than about 4 m are parameterized using a modified Smagorinsky model (Kaltenbach et al., 1994). Subgrid-scale tracer fluxes also depend on a variable subgrid-scale Prandtl number, which is parameterized in terms of the grid-scale gradient Richardson number as described in Whitt and Taylor (2017). Subgrid-scale Schmidt numbers associated with biogeochemical tracers are set equal to the subgrid-scale Prandtl number.

Two configurations of the BLES are compared, a frontal zone configuration (e.g., Callies & Ferrari, 2018;Skyllingstad et al., 2017;Taylor & Ferrari, 2010;Thomas & Taylor, 2010; LES,F) with submesoscales and a configuration without a front (LES,NF) where submesoscales are excluded. The LES,NF configuration captures the conventional effects of small-scale turbulence (e.g., McWilliams & Sullivan, 2000;Skyllingstad et al., 2000;Pham & Sarkar, 2017). Differences between the two configurations quantify how submesoscales modify these conventional effects. In the LES,F configuration, the flow is expressed as a periodic (in x and y) perturbation from a fixed mean horizontal density gradient

and thermal wind shear ⟨M 2 ⟩ x,y ∕f = 5 × 10 -4 s -1 that are qualitatively representative of the OSMOSIS site during September 2012. Here, ⟨⟩ x,y denotes a horizontal average, 𝜌 is the density, 𝜌 0 = 1,026 kg/m 3 is the reference density, g is the acceleration due to gravity, the buoyancy b =-g𝜌∕𝜌 0 , and f = 10 -4 s -1 is the Coriolis frequency. The chosen ⟨M 2 ⟩ x,y = 5 × 10 -8 s -2 is representative of roughly the 75th percentile of all horizontal buoyancy gradients |∇ h b| observed at 2 to 5 km length scales by ocean gliders during autumn at the OSMOSIS site (Thompson et al., 2016). This means the chosen horizontal density gradient is strong but typical of the region and season; stronger gradients are expected roughly a quarter of the time and weaker gradients are expected roughly three quarters of the time. Hence, although we sometimes refer to the mean horizontal density gradient as a "front," the gradient is thought to be fairly typical of the open ocean and, in particular, is 2 orders of magnitude weaker than strong submesoscale fronts that are occasionally observed in boundary currents (D'Asaro et al., 2011).

The turbulent physical initial condition at the onset of the storm with finite amplitude submesoscale eddies (Figure 1) is obtained from a 3-day spin-up simulation that is initialized with low-amplitude red noise in the frontal zone with a vertical density profile based on the observed density profile above the Porcupine Abyssal Plain during September 2012. The spin-up simulation is forced by a constant air-sea buoyancy flux B A = 3 × 10 -9 m 2 /s 3 (equivalent to about 10 W/m 2 cooling) that is chosen to balance the restratifying effects of submesoscale mixed layer baroclinic instabilities, as discussed in Whitt and Taylor (2017; see also Mahadevan et al., 2010Mahadevan et al., , 2012)). At the start of the spin-up simulation, the fastest-growing mixed layer baroclinic instability has a horizontal length scale of about 985 m and a growth timescale of about 13 hr, as discussed in Whitt and Taylor (2017; see also Stone, 1966). The domain size is chosen so that it contains two wavelengths of the fastest growing baroclinic mixed layer instability mode in each horizontal dimension. Hence, the domain is 1,970 m by 1,970 m by 80 m and covered by a grid with 1,024 by 1,024 by 160 points that achieves a uniform resolution of 1.9 m by 1.9 m by 0.5 m in x (along-front) and y (across-front) and z (vertical), respectively. As a result, the initial perturbations grow and the physical fields have strong submesoscale variability at the onset of the storm, as shown in Figure 1.

The LES,NF configuration is identical to the LES,F configuration except that the domain is smaller (492.5 m by 492.5 m by 80 m) and ⟨M 2 ⟩ x,y = 0. In particular, this analog configuration uses the same code, the same grid resolution, the same surface boundary conditions, and the same mean density profile ⟨𝜌⟩ x,y (z) at day 0 as the large domain. But there is no initial submesoscale variability. Rather, the LES,NF configuration is initialized with low-amplitude white noise in the velocity components and no initial density perturbations.

Storm Forcing

The storm forcing is represented by an idealized spatially uniform but time-dependent surface stress that points 45 • to the right of the geostrophic shear; it is inspired by an observed storm during 24-26 September 2012 at the OSMOSIS site (Rumyantseva et al., 2015). The idealized stress ramps up linearly over 2 days to a maximum amplitude of 0.6 N/m 2 and then ramps down linearly to zero during the third day (as shown in Figure 2). Sensitivity studies with a maximum stress of 0.3 N/m 2 are qualitatively similar and presented in the supplemental material (see Table 1 for a full list of simulations). Unless otherwise noted, the text refers to the cases with the stronger winds. The direction, magnitude, and duration of the wind forcing are chosen to be similar to the observed amplitude, direction, and duration of the wind stress during the storm based on NCEP reanalysis winds, satellite sea surface temperature imagery, and results reported in Rumyantseva et al. (2015) and Buckingham et al. (2017).

In the presence of a mean horizontal buoyancy gradient at the surface (as in LES,F), the Ekman flow can produce a buoyancy tendency in the mixed layer via horizontal advection. For example, if the wind stress includes a down-front component, meaning a component in the direction of the thermal wind shear k∕f ×∇ h b, then the Ekman flow tends to reduce the mixed layer buoyancy, energize submesoscale symmetric instabilities aligned parallel to the front, and may cause entrainment. But, if a component of the wind is up front, meaning opposite to the direction of the thermal wind shear, then the Ekman flow tends to increase the mixed layer buoyancy and may cause subduction (Brannigan, 2016;Whitt et al., 2017;Whitt et al., 2017). In the simulations presented here, the wind includes a down-front component, and the horizontally averaged Ekman buoyancy flux EBF= 𝜏 x ⟨M 2 ⟩ x,y ∕𝜌 0 f, which is driven by the down-front component of the wind stress 𝜏 x , reaches maximum magnitudes of 2 × 10 -7 m 2 /s 3 (or about 660 W/m 2 Ekman heat flux) when the magnitude of the stress is 0.6 N/m 2 (Figure 2). It may be noted that the horizontally averaged EBF is independent of the interior evolution of the simulations, because the stress is a spatially uniform surface boundary condition (and therefore uncorrelated with any interior spatial variability of density), and the mean horizontal density gradient ⟨M 2 ⟩ x,y is an externally imposed constant.

Following the storm, the simulations continue for about 3 days without wind stress. For simplicity, the air-sea buoyancy flux is held constant at B A = 3 × 10 -9 m 2 /s 3 (or about 10 W/m 2 heat loss) during and after the simulated storm.

It may be noted that the simulations presented here do not have realistic air-sea buoyancy loss, which is usually enhanced during storms. That is because, although we have run simulations with more realistic buoyancy loss (not shown), the stronger buoyancy loss has little impact on the results. A wind stress of 0.6 N/m 2 is a far stronger source of energy for turbulent entrainment in a 40-m-deep mixed layer than the roughly 200 W/m 2 net heat loss observed by Rumyantseva et al. (2015), as can be inferred from scaling arguments (Kraus & Turner, 1967;Pollard et al., 1972) and observations (Thomas et al., 2016). Further, even a realistic surface buoyancy flux would be about a factor of 3 smaller than the EBF. Using the observationally supported scalings reported in Thomas et al. (2016), the depth-integrated dissipation of turbulent kinetic energy should have contributions of about 6 × 10 -5 m 3 /s 3 from wind and 8 × 10 -6 m 3 /s 3 from EBF during peak wind. That is to say, the direct injection of energy by the wind during the storm is the dominant source of energy for the turbulence. Based on these scalings, the turbulence dynamics are not expected to differ much between the various scenarios presented here, nor is the biogeochemical response.

Biogeochemical Model

The biogeochemical model is a modified version of the NPZD (nutrient, phytoplankton, zooplankton, and detritus) model implemented by Whitt, Taylor, et al. (2017) and Whitt, Lévy, et al. (2017), which is based on classic models like that of Fasham et al. (1990) and a previous NPZD implementation by Powell et al. (2006).

The biogeochemical model equations are

where D∕Dt = 𝜕∕𝜕t + u •∇is the material derivative, the subgrid-scale diffusivity is given by 𝜅 SGS (which varies in time and space), the molecular diffusivity is set to a constant 𝜅 0 = 10 -6 m 2 /s, and the constant parameters are listed in Table 2.

The biogeochemical initial conditions (Figure 3), which are imposed as horizontally uniform fields at t = 0 d, are from an equilibrium (i.e., steady state) solution to the biogeochemical model equations with the parameters in Table 2 and a constant vertical diffusivity 𝜅 SGS +𝜅 0 = 2×10 -5 m 2 /s. The initial Z and D profiles (not shown) have a similar vertical structure to P, with maximum concentrations of 0.22 and 0.35 mmol/m 3 , respectively. This equilibrium solution is calculated numerically by running one-dimensional simulations of the biogeochemical model equations for 10 years with a one-dimensional NPZD model described earlier (Whitt, 2017;Whitt, Taylor, et al., 2017;Whitt, Lévy, et al., 2017). In addition, numerical experimentation suggests that the resulting equilibrium profile is effectively independent of the initial conditions of N, P, Z, and D for the parameter set in Table 2.

Journal of Geophysical Research: Oceans 10.1029/2019JC015370

Table 2 Biogeochemical Model Parameters Parameter Symbol Value Light attenuation due to seawater k w 0.055 m -1 Initial slope of the productivity-irradiance curve 𝛼 0.1 d -1 (W/m 2 ) -1 Surface photosynthetically available radiation I 0 50 W/m 2 Phytoplankton max. uptake rate V m 1.5 d -1 Half-saturation for nutrient uptake by phy. k N 0.8 mmol/m 3 Phytoplankton mortality rate to detritus 𝜎 d 0.04 d -1 Ivlev constant for zooplankton grazing 𝛬 0.5 (mmol/m 3 ) -1 Maximum zooplankton grazing rate R 0.9 d -1 Zooplankton excretion to nutrient 𝛾 N 0.3 Linear zooplankton mortality rate 𝜁 0.04 d -1 Quadratic zooplankton mortality rate ζ 0.28 d -1 (mmol/m 3 ) -1 Detritus remineralization rate 𝛿 0.1 d -1 Detrital sinking velocity w d -0.25 m/d Restoring nutrient concentration at 80 m N 80 2.5 mmol/m 3 Inverse restoring timescale for nutrient at 80 m 𝛽 N 0.25 d -1

Due to the construction of the biogeochemical initial condition, the horizontally averaged nutrient ⟨N⟩ x,y and phytoplankton ⟨P⟩ x,y profiles have strong vertical gradients throughout the top 60 m at Day 0 ( Figure 3).

In the upper 30 m, these strong vertical biogeochemical gradients contrast with weak vertical density gradients, which are relatively permissive of turbulent mixing compared to the pycnocline below. Hence, shortly after the simulations begin, storm-driven turbulence drives rapid vertical mixing and homogenization of the mean biogeochemical concentrations over the top 30 m.

To reflect the idealized nature of the simulations, we refer to the simulated nutrient as "nutrient" throughout the paper. In addition, it is not our intention to validate the model or exactly reproduce observations here.

However, the biogeochemical model parameters in Table 2 are tuned within reasonable ranges (see the appendix) so that the chosen initial equilibrium, which reflects the parameters, has a nutrient profile that is similar to the observed dissolved silica profile in autumn 2012 above the Porcupine Abyssal Plain on Journal of Geophysical Research: Oceans 10.1029/2019JC015370

21-22 September 2012 a few days before the storm (Figure 3b). The subsurface maximum in the initial phytoplankton profile is also qualitatively consistent with observations of vertical profiles of chlorophyll fluorescence (Figure 3a).

Column Setup in ROMS

In order to assess the effectiveness of some existing ocean model parameterizations at representing the horizontally averaged behavior of LES, results from the BLES models are compared to results from a one-dimensional column model implemented in the regional ocean modeling system (ROMS+NPZD) (e.g., Powell et al., 2006;Whitt, Taylor, et al., 2017). Two ROMS configurations are used that can be viewed as parameterized analogs of the two LES configurations, one with parameterized submesoscales (ROMS,F) and one with no submesoscales (ROMS,NF). In the ROMS,NF configuration, the physical effects of turbulence are parametrized using the conventional KPP mixing scheme (Large et al., 1994), as implemented in the public version of ROMS (myroms.org). There are no effects of submesoscales and no horizontal density gradient, that is, ⟨M 2 ⟩ x,y = 0. In the ROMS,F configuration, the physical effects of both turbulence and submesoscales are parametrized rather than explicitly resolved, and the same mean horizontal density gradient as in the LES,F experiments is used, that is,

In all column model experiments, ROMS solves for the time evolution of horizontal momentum, temperature (which is equivalent to density via a linear equation of state), and biogeochemical tracers on the same uniform vertical grid with 160 levels and 0.5 m resolution that is used in the LES. Vertical fluxes are computed using a conservative parabolic spline method. In addition, the column model is forced by the same surface momentum and buoyancy fluxes and solves the same biogeochemical model equations.

The parameterizations for small-scale turbulence and for submesoscales are implemented in the column version of ROMS as follows. First, the KPP surface layer scheme is modified so that the boundary layer depth and thereby the diffusivity and viscosity profiles are defined based on a bulk Richardson number that depends on the constant geostrophic shear ⟨M 2 ⟩ x,y ∕f as well as the evolving ageostrophic shear. Second, terms are added to the right side of the temperature equation (which is also the density equation since we are using linear equation of state) to represent horizontal advection of the background lateral temperature gradient, including an extra term of the form v⟨M 2 ⟩ x,y ∕g𝛼, where 𝛼 is a constant thermal expansion coefficient. This extra advection term includes the EBF defined earlier. In addition, a second term is added to the right-hand side of the temperature equation of the form v * ⟨M 2 ⟩ x,y ∕g𝛼, where v * is defined below. This second new term in the temperature equation represents the restratifying effect of submesoscale mixed layer baroclinic instabilities via the parameterization of Fox-Kemper et al. (2008;FFH), in which the temperature tendency associated with this restratifying effect is expressed as -∇ • (u * T), where

and the other components of the bolus velocity vector u * are equal to zero in this case. The time-dependent mixed layer depth MLD 3 in (2) is defined in section 2.2.2 below. The time-dependent vertical structure function 𝜇(z), which depends on MLD 3 as well as z, is given by (39) of FFH. The coefficient C = 0.06, which FFH find to be the best fit in a series of simulations without diurnal convection or wind, but this parameter C is subject to some uncertainty (e.g., Callies & Ferrari, 2018;Bachman & Taylor, 2016;Taylor, 2016). Finally, the convergence of the viscous vertical flux of geostrophic momentum is added to the tendency. This flux of geostrophic momentum is calculated like the viscous flux of ageostrophic momentum, by multiplying the same vertically and temporally variable viscosity coefficients to the constant background geostrophic shear profile.

Model Diagnostics 2.2.1. Submesoscales and Smaller-Scale Turbulence: Scale Separation

Throughout this paper, the term submesoscale is used descriptively to refer to features of the velocity, density, and tracer fields that reside primarily in the ocean mixed layer with characteristic horizontal length scales of 0.1-10 km and, unless stated otherwise, does not imply any particular physical process. In the context considered in this paper-the midlatitude North Atlantic during autumn-this definition implies that submesoscales have horizontal scales that are larger than the mixed layer depth, which is generally 10-100 m. Therefore, submesoscales have small aspect ratios and are distinct from conventional mixed layer turbulence, which is typically composed of features with horizontal wavelengths of the same order as or smaller

Journal of Geophysical Research: Oceans 10.1029/2019JC015370

than the mixed layer depth (i.e., ≲0.1 km horizontal scales) and order-one aspect ratios. Since some simulations reveal a change in the slope of kinetic energy and tracer spectra at a wavelength slightly larger than the mixed layer depth (Hamlington et al., 2014;Smith et al., 2016;Whitt & Taylor, 2017), we define a separation at the horizontal wavelength of 150 m to separate smaller turbulence from larger submesoscales in the computation of the variances and vertical fluxes (the rationale behind this particular choice is discussed in Whitt & Taylor, 2017). For example, the total vertical nutrient flux ∫ ŵ(k) N * (k) dk is separated into its submesoscale and turbulent parts, using a separating wavenumber 1∕150 cyc/m. The results presented here do not change qualitatively if this separation wave number is increased or decreased by a factor of 2. Hence, we say the LES,F scenarios in the larger 2 km domain contain submesoscales, and we say the LES,NF scenarios in the smaller 0.5 km domain do not contain submesoscales (see Table 1 for a list of simulations).

Potential Energy Change Due to Vertical Mixing

The potential energy (PE = ∫ 𝜌g(z + H) dz where z ranges from -80 mt o0ma n dH = 80 m is the domain depth (e.g., Burchard & Hofmeister, 2008;Crawford & Large, 1996) can be increased in three ways in response to the storm. First, vertical mixing acts to raise the center of mass and thereby increases the potential energy. Second, vertically sheared horizontal mean currents can act on the mean horizontal density gradient to reduce the vertical stratification and hence also increase the potential energy. For wind-driven (Ekman) currents, the corresponding buoyancy flux is referred to as an EBF. Third, PE can increase due to a

Journal of Geophysical Research: Oceans 10.1029/2019JC015370

stabilizing air-sea buoyancy flux (e.g., surface heating). See Whitt and Taylor (2017) for more details on the energetics.

In order to quantify the approximate increase in PE due to mixing, the change in PE due to lateral advection and the air-sea buoyancy flux, B A , must be subtracted from the change in PE, ΔPE, that is,

(4)

Dividing the residual change in potential energy by the wind work ΔPE * ∕ΔKE, where ΔKE = ∫ ⟨•u dt and u is the surface ocean velocity, approximately quantifies the fraction of the wind work that goes toward increasing the potential energy via turbulent mixing, that is, a mixing ratio r m =ΔPE * ∕ΔKE.

Density

At the onset of the storm, density is already fairly well mixed in the upper 30 m (Figure 3). Yet the storm is sufficiently strong to mix denser water up from the pycnocline, to deepen MLD 3 by approximately 10 m and to increase the surface density (Figures 4c and 4f). The deepening of MLD 3 is approximately the same in the scenarios with and without submesoscales, but the MLD 3 deepens roughly half as much when the wind stress is half as strong (Figure S5). That is, the deepening of MLD 3 is insensitive to the horizontal density gradient or the presence of submesoscale heterogeneity, up to the degree each are present in these simulations. However, this is not the case for MLD 05 , which shoals during the storm with submesoscales and deepens without submesoscales (Figure 4; see also Figure S5). The MLD 05 is shallower in the simulation with submesoscales because the vertical density gradient above MLD 3 is sensitive to the presence of the submesoscales during the storm. Without submesoscales, the vertical density gradient above MLD 3 is concentrated very close to MLD 3 , whereas with submesoscales the vertical density gradient above MLD 3 is distributed more evenly with depth. Hence, the curvature of the vertical density profile is weaker near MLD 3 with submesoscales, which is consistent with enhanced vertical mixing of density across MLD 3 during the storm in that case. We can also note that the surface density is greater after the storm with submesoscales but this denser surface water could be attributable to the horizontal EBF and/or enhanced vertical transport from the pycnocline. Therefore, the denser surface water with submesoscales is not necessarily associated with increased vertical mixing.

In order to show that the submesoscale heterogeneity triggers enhanced vertical exchange of density from the pycnocline to the mixed layer during the storm, we quantify the increase in potential energy due to vertical mixing alone (ΔPE * , defined in section 2.2.5). We find that ΔPE * is about twice as large with submesoscales as without ( Figure 5), and we obtain a quantitatively similar result for both wind stress magnitudes (Figure S5). The mixing ratio (r m =Δ PE * ∕ΔKE, defined in section 2.2.5) ranges from about 1.5% without submesoscales to 3% with submesoscales (Figure 5b), and r m is only reduced by about 1/4 when the stress is half as strong, although both ΔPE * and the kinetic energy input by the wind ΔKE are reduced by about 2/3 with the weaker stress (Figure S6). These results show that the submesoscales enhance the net (sub-mesoscale+turbulent) vertical transport of density across the mixed layer base relative to the simulations without submesoscales, without significantly increasing the kinetic energy input by the wind.

Before discussing the biogeochemical response, it is worth noting that a parameterization of r m is a crucial explicit or implicit choice (e.g., via a bulk Richardson number threshold, as discussed in Pollard et al., 1972) in many models of ocean boundary layer mixing and stratification (Burchard & Hofmeister, 2008;Crawford & Large, 1996;Kraus & Turner, 1967;Large et al., 1994;Pollard et al., 1972), and these choices can strongly impact the modeled response of the mixed layer to wind forcing. However, this is the first time that LES Journal of Geophysical Research: Oceans 10.1029/2019JC015370 Figure 6. Full-depth column integrals of phytoplankton ⟨P⟩ x,y (a) and storm-driven ⟨NCP⟩ x,y (b), averaged over three time periods. The reduction in column-integrated ⟨N⟩ x,y is equal to the storm-driven ⟨NCP⟩ x,y here.

results have been used to explicitly show that submesoscale heterogeneity can systematically enhance r m , and no mixing scheme that is used operationally in ocean models includes a parameter to explicitly enhance the mixing efficiency to account for unresolved submesoscale effects like those simulated here. Based on this new result, it is to be expected that submesoscales also enhance the mean vertical nutrient flux ⟨wN⟩ x,y and net community production ⟨NCP⟩ x,y in this nutrient-limited autumn storm scenario. In what follows, we quantify and explore these biogeochemical impacts in more detail.

Nutrient and Phytoplankton

We start by examining how the horizontally averaged phytoplankton and nutrient profiles evolve during and after the storm in the LES,F and LES,NF simulations. By some measures, the mean biogeochemical response to the storm is similar with and without submesoscales (Figure 4). Storm-driven turbulence fluxes nutrients and phytoplankton up toward higher light levels near the ocean surface. The storm-driven upward biogeochemical fluxes trigger a reduction in column-integrated nutrient ⟨N⟩ x,y and a nearly equal and opposite increase in column-integrated phytoplankton biomass ⟨P⟩ x,y (Figure 6). This mirroring of the changes in column-integrated ⟨N⟩ x,y and ⟨P⟩ x,y indicates that primary productivity is the dominant mechanism of column-integrated biogeochemical transformation. These results are qualitatively similar under the weaker stress, but the changes in column-integrated ⟨N⟩ x,y and ⟨P⟩ x,y are about 1/3 to 1/2 as large, which is qualitatively consistent with less deepening of the MLD 3 and less nutrient entrainment under weaker wind stress (Figure S5). In addition, this evolution is qualitatively consistent with expectations based on available observations of the nutrient and/or chlorophyll responses to autumn storm passage (e.g., Babin et al., 2004;Rumyantseva et al., 2015).

However, the mean biogeochemical response to the storm also reveals some systematic differences with and without submesoscales. First, submesoscales enhance the column-integrated ⟨NCP⟩ x,y by about a factor of 2, which is qualitatively consistent with the increased mixing ratio r m . In addition, the column-integrated ⟨NCP⟩ x,y is also enhanced by submesoscales when the wind is half as strong, albeit by somewhat less than a factor of two (Figure S5). We can also note a reduced curvature of the mean vertical nutrient and phytoplankton profiles near MLD 3 in the scenario with submesoscales, which is consistent with the reduced curvature of the density profiles near MLD 3 and enhanced vertical mixing there with submesoscales. Finally, we also find faster surface accumulation of phytoplankton and stronger vertical gradients in both phytoplankton and nutrient that develop within the mixed layer in the scenarios with submesoscales compared to the scenarios without submesoscales. In the following section, we show how submesoscales modulate the net (submesoscale+turbulent) vertical fluxes of biogeochemical constituents to explain these results.

Horizontal Variability in the LES Solutions

Based on the results of the previous section, it is expected that average vertical nutrient fluxes ⟨wN⟩ x,y will be sensitive to the presence of submesoscales. However, since the biogeochemical dynamics, such as primary productivity and grazing, are nonlinear, the biogeochemical covariability at scales from meters to submesoscales, such as cov(N, P), may rectify on larger-scale biogeochemical dynamics as well. In addition, the fluxes and covariances are best understood in the context of the variances, such as var(N), on which they depend. Hence, this section addresses the following questions for the first time with BLES: Does submesoscale heterogeneity impact the horizontal biogeochemical variances, such as var(N), or key covariances, such as cov(N, P), which impacts primary productivity, or cov(w, N)=⟨wN⟩ x,y , in this storm scenario? To what degree do the submesoscales (>150 m) and turbulent scales (<150 m) contribute to the total variances and covariances during the different phases of the storm? Is there a rectified effect of submesoscales on physical and biogeochemical variances and covariances at turbulent scales?

We find that the physical/biogeochemical variances and covariances are greatly modified by the presence of submesoscales during and after the storm. To gain further perspective on this aspect of the solutions, we first qualitatively compare the lateral variability of the nutrient N and nutrient flux wN in snapshots from the different solutions during the storm (Figure 7) before proceeding to statistical measures of the variability. Just above MLD 3 (i.e., MLD 3 + 4 m), the variability in nutrient concentration differs qualitatively between LES,F and LES,NF. With submesoscales, the nutrient distribution reflects the submesoscale vertical displacements of the nutricline (Figures 7a), which dominate the total variability. However, the dominant characteristic scales of the nutrient flux wN (which are qualitatively similar to w, not shown) are much smaller. Yet these dominant small-scale structures in wN are modulated somewhat by the larger submesoscales ( Figure 7b). In contrast, the LES,NF scenario exhibits no obvious qualitative evidence of organization at submesoscales in either N or wN (Figure 7), and the characteristic horizontal scales of N and wN variability are comparable in magnitude to MLD 3 . In this sense, the simulated variability of N and wN in the LES,NF scenario is qualitatively consistent with classic ideas about boundary layer turbulence, in which the largest aspect-ratio-one turbulent eddies with a length scale similar to MLD 3 dominate both the variances and the fluxes. In contrast, the LES,F scenario is qualitatively inconsistent with classic ideas about boundary layer turbulence in the same sense, because low-aspect-ratio structures that are a factor of 10-50 times larger than MLD 3 qualitatively dominate the variance (for a spectral perspective on this statement, compare the schematic power spectrum of boundary layer turbulence in Figure 1 of Wyngaard, 2004, in which the power drops at horizontal wavelengths greater than the boundary layer height, to simulated power spectra in Figure 5 of Whitt & Taylor, 2017, in which the power rises at horizontal wavelengths greater than the mixed layer depth).

Biogeochemical Variances in the LES solutions

For a more statistical perspective on the dependence of the biogeochemical variability on the submesoscales, we explicitly quantify and compare the vertical profiles of the nutrient variance var(N), and we show how submesoscales (>150 m) and smaller turbulent scales (<150 m) contribute to the total variance in the LES,F and LES,NF scenarios (Figure 8). The phytoplankton variance profiles var(P) are qualitatively similar to var(N), and the results are qualitatively insensitive to reducing the wind stress by half (see Figures S7 and S8), so we only discuss the stronger wind scenario in this context.

During the storm (Days 2 and 3), var(N) is maximal just below MLD 3 and declines monotonically upward so that it is 3 to 4 orders of magnitude smaller at the surface. With submesoscales, there is much greater total variance in the mixed layer due almost entirely to submesoscale (>150 m) variance, similar to weakly forced LES simulations of passive tracers released in the pycnocline by Smith et al. (2016). However, in a more novel Journal of Geophysical Research: Oceans 10.1029/2019JC015370 result, the presence of submesoscales also increases the variance of N at smaller (<150 m) turbulent scales during the storm by up to about an order of magnitude, particularly in the upper part of the mixed layer (Figure 8c). The enhanced smaller-scale N variance with submesoscales is likely due primarily to enhanced physical variance production from stronger mean and submesoscale vertical gradients (Figure 4), since the smaller-scale turbulent kinetic energy and the dissipation of turbulent kinetic energy in the mixed layer are only modestly modified by submesoscales during the storm (see Figure S12 here and Figure 5 of Whitt & Taylor, 2017).

Although kinetic energy at submesoscales and turbulent scales declines substantially after the storm (Whitt & Taylor, 2017; see also Figure S12), var(N) generally increases after the storm and remains about 1 to 2 orders of magnitude greater with submesoscales than without. With submesoscales, var(N) achieves its maximum in the of the mixed layer, that is, from about 10-to 30-m depth (Figure 8f), due to a substantial increase in submesoscale variance at these depths after the storm. This increase in submesoscale variance after the storm is likely due primarily to the near shutdown of turbulent mixing and the corresponding reduction of the rate of dissipation of submesoscale N variance via turbulent mixing. On the other hand, production of submesoscale variance is not substantially enhanced after the storm, since submesoscale kinetic energy declines after the storm (Whitt & Taylor, 2017), and mean vertical tracer gradients are only slightly enhanced with submesoscales (Figure 4).

Nutrient Fluxes in the LES Solutions

Here, we quantify and discuss the vertical nutrient flux, and its submesoscale and turbulent parts in the LES,F and LES,NF scenarios (Figure 9; var(N) is in Figure 8; var(w) is in Figure S12; and multiscale corrrelation profiles corr(w, N) are in Figure S13). Consistent with the results above, we find that both the total fluxes and their respective submesoscale and turbulent parts differ by an order-one fraction with and without submesoscales, both during and after the storm.

First, we find that the presence of submesoscales substantially enhances the total vertical nutrient flux during the storm (between Days 2 and 3, left column in Figure 9). In addition, the effect of the submesoscales on the maximum total flux is qualitatively similar but greater in the scenario with stronger stress (a 115% increase) than the scenario with the weaker stress (a 34% increase; see Figure S10). However, the total flux profiles exhibit a qualitatively similar vertical structure in scenarios with and without submesoscales. In both scenarios, the maximum flux occurs between about 0.8MLD 3 and 0.9MLD 3 , and the flux drops almost linearly to zero between the depth where it achieves its maximum and the surface above or about 1.1MLD 3 below. Yet the strong fluxes penetrate deeper below the mixed layer MLD 3 in the LES,F scenarios compared to the LES,NF scenarios, which is consistent with deeper penetration of high rates of dissipation of turbulent kinetic energy (Whitt & Taylor, 2017) and vertical velocity variance (Figure S12) into the pycnocline/nutricline in the presence of submesoscales.

The partitioning of the flux between submesoscales (>150 m) and smaller turbulent scales (<150 m) during the storm differs significantly depending on the presence or absence of the submesoscales, and the differences are qualitatively similar with weaker wind (Figure S10). Without submesoscales in LES,NF, the total flux is dominated by smaller turbulent scales. On the other hand, with submesoscales in LES,F, submesoscale and smaller-scale turbulent fluxes are both positive and submesoscales contribute an order-one fraction of the total flux. The strong submesoscale flux is qualitatively consistent with strong submesoscale N variance (Figure 8), and it shows that the submesoscale heaving of the nutricline observed in Figure 7 is also associated with some irreversible upward nutrient transport at submesoscales. With a weaker wind, the fraction of the total flux from submesoscales is larger, but the magnitude of the submesoscale flux is smaller (see Figure S10). Incidentally, the decrease in submesoscale flux with the decrease in wind stress magnitude is qualitatively consistent with previous process studies with ocean models that parameterize small turbulent scales (e.g., Capet et al., 2008;Mahadevan & Tandon, 2006;Thomas et al., 2008;Whitt, Taylor, et al., 2017), although the causes of there relationships are not necessarily the same (see Whitt & Taylor, 2017, for some relevant discussion). However, in a more novel result, the smaller-scale turbulent contribution to the nutrient flux is significantly greater during the storm with submesoscales than without, even though vertical velocity variance is only modestly altered throughout most of the mixed layer (Figure S12). This enhanced flux (i.e., covariance) is qualitatively consistent with enhanced turbulent N variance with submesoscales, but the correlation between nutrient and vertical velocity is actually smaller, both in aggregate and at turbulent scales, in the simulations with submesoscales than without submesoscales (Figure S13). Finally, perhaps the most intriguing feature of the nutrient flux profiles during the storm is the deeper penetration of the strong vertical fluxes below the mixed layer in the LES,F scenarios with submesoscales compared to the LES,NF scenarios without. This enhanced flux near MLD 3 occurs in conjunction with a substantial increase in smaller-scale turbulent vertical velocity variance (Figure S12) and is due mostly although not entirely to stronger smaller-scale fluxes (Figure 9). After the storm, the nutrient flux profiles change dramatically compared to their counterparts during the storm. In particular, the magnitude of the total flux is reduced by about an order of magnitude in both simulations (note the different x axes in the left and right columns in Figure 9). In addition, the maximum total flux is about a factor of 2 weaker with submesoscales than without. The partitioning of the total flux between submesoscales and turbulent scales is also very different between the two scenarios. In particular, the smaller-scale turbulent flux is much weaker in the scenario with the submesoscales (LES,F) than in the scenario without submesoscales (LES,NF). Consequently, after the storm, the total flux is achieved mostly by the submesoscale flux in LES,F and by the smaller-scale turbulent flux in LES,NF. The shutdown of the turbulent flux after the storm in the LES,F scenario and increase in vertical biogeochemical gradients is qualitatively consistent with earlier results showing how submesoscales influence the timing of the spring bloom by enhancing the density stratification, suppressing turbulence, and thereby modifying vertical fluxes in deeper spring mixed layers (Taylor, 2016). More generally, these results suggest that fronts/submesoscales may contribute to the rapid formation of biogeochemical gradients within mixed layers after strong wind events during any time of year by suppressing turbulent fluxes even when the mixed layer depths have not yet shoaled significantly (e.g., as observed by Carranza et al., 2018).

Scenarios With a Front and Submesoscales

The frontal zone version of the column model ROMS,F does not perform as well at capturing the physical/biogeochemical dynamics of the LES,F scenario with the front/submesoscales. In fact, the ROMS,F scenarios are in many ways more similar to the ROMS,NF and LES,NF scenarios than their analog LES,F scenarios, which shows that the submesoscale heterogeneity simulated by LES,F is the cause of the differences relative to LES,NF, not the front. As in previous sections, only the scenarios forced by the stronger winds are presented here; the scenarios forced by the weaker winds are qualitatively similar and appear in Figures S2-S6 and S10-S11.

The most distinctive result emerging from the LES,F/LES,NF comparison above is that the presence of submesoscales in LES,F roughly doubles the vertical nutrient fluxes and mixing ratio during the storm, and the submesoscales indirectly cause turbulence and turbulent fluxes to penetrate below MLD 3 . In contrast, the total vertical nutrient flux in the ROMS,F scenario is actually reduced compared to the ROMS,NF and LES,NF scenarios during the storm. As a result, the maximum vertical nutrient flux in ROMS,F is about a quarter of that in LES,F (Figure 10a). Further, ROMS,F fails to capture the deeper penetration of vertical fluxes below MLD 3 that is simulated in LES,F.

Comparing the vertical diffusivity profiles 𝜅 N during the storm (Figure 10c) shows that ROMS,F fails to capture the strong vertical nutrient fluxes in the mixed layer partly because the mean nutrient gradients are too weak, not just because the diffusivities are too weak. For example, between about 0.5MLD 3 and 0.9MLD 3 , the diffusivity is actually lower in LES,F than ROMS,F, while the fluxes are about three to four times larger in LES,F than ROMS,F (mean nutrient profiles are compared in Figure S3). On the other hand, just below the mixed layer (between about MLD 3 and 1.1MLD 3 ), both the flux and diffusivity are much larger in the LES,F than ROMS,F. These results suggest that submesoscales systematically modify the empirical nondimensional shape functions as well as the turbulent boundary layer depth, which characterize the diffusivity profiles in the mixing scheme (Large et al., 1994), during the storm.

After the storm, there is some indication that the physics encapsulated in the ROMS,F scenario may capture some of the essential physics in the LES,F scenario. In particular, the maximum total fluxes in the mixed layer are smaller in magnitude and occur at a shallower depth than in the corresponding LES,NF and ROMS,NF scenarios (cf. Figures 10b and 11b). In addition, both the mixed layer depth MLD 05 and the surface mixing layer depth shoal in the ROMS,F scenario. However, there are still some discrepancies between LES,F and ROMS,F. For example, the maximum nutrient flux is only about a quarter to a third of that in LES,F. In the upper half of the mixed layer, where vertical diffusivities in ROMS,F are mostly larger than in LES,F, the reduced vertical fluxes are caused by too-weak mean gradients. However, between 0.5MLD 3 and MLD 3 , the reduced vertical fluxes are mostly caused by too-small vertical diffusivities. In addition, the reduced fluxes are presumably partially a result of missing submesoscale vertical fluxes, which dominate in LES,F after the storm, but are not included in ROMS,F. Despite differences, these significantly reduced diffusivities after the storm-in both ROMS,F and LES,F-manifest in enhanced mean vertical gradients in phytoplankton above MLD 3 by Days 5 and 6 in the scenarios with a front compared to their counterparts without a front, that is, LES,NF and ROMS,NF (Figures S2 and S3). Hence, the relatively rapid formation of biogeochemical vertical gradients within the mixed layer in both LES,F and ROMS,F scenarios is partly attributable to the substantial reduction in vertical mixing in the scenarios in a frontal zone. However, as a result of the failure of the ROMS,F column model to capture the enhanced vertical nutrient fluxes during the storm, there is only a third to a half as much net community production and accumulation of column-integrated phytoplankton during and shortly after the storm (Figure 12b).

Summary of the Comparisons Between the Column Model and LES

Qualitatively and quantitatively, the conventional boundary layer column model ROMS,NF accurately captures the physical-biogeochemical evolution of the conventional boundary layer LES,NF without a front/submesoscales. This is encouraging, with the caveat that there is some indication that the performance of KPP may not be equally good in other parts of physical parameter space (Pham & Sarkar, 2017). Regardless, the column model with a frontal zone ROMS,F fails to reproduce the physical-biogeochemical evolution in the LES with a frontal zone. In particular, the column model fails to capture the fact that submesoscale physical variability in the mixed layer enhances vertical nutrient fluxes, phytoplankton accumulation, and Journal of Geophysical Research: Oceans 10.1029/2019JC015370 net community production during the life cycle of the storm. However, the column model does capture the qualitatively reduced mixing within the mixed layer MLD 3 and the related rapid emergence of vertical phytoplankton gradients within the mixed layer after the storm with a front/submesoscales compared to scenarios without a front/submesoscales.

The failure of the column model to capture the enhanced vertical mixing and nutrient fluxes during the storm rules out several a priori plausible explanations for this difference between LES,F and LES,NF, including the EBF, vertical mixing of geostrophic momentum, and/or the restratifying buoyancy flux associated with mixed layer baroclinic instabilities as described by FFH via (2). However, these three factors could contribute indirectly, in combination with other properties of submesoscale variability in LES,F, to induce enhanced mixing during the storm. On the other hand, the partial success of the column model in capturing the reduced vertical mixing and the more rapid formation of vertical gradients in phytoplankton after the storm suggests that the column model captures some of the essential physical-biogeochemical dynamics associated with the rapid formation of these vertical phytoplankton gradients in LES,F. In particular, the combination of the EBF, mixing of geostrophic momentum, and restratification via (2) results in substantially reduced mixing within the mixed layer after the storm, and the more rapid formation of biogeochemical gradients within the mixed layer. Future work is required to generalize these results and to develop and test other hypotheses about the key physical-biogeochemical mechanisms responsible for the enhanced vertical mixing and net community production during the storm.

Summary

This study uses process-oriented numerical simulations to address the questions of whether and to what degree submesoscale heterogeneity influences the physical-biogeochemical response to a storm. Previous process studies using numerical models with parameterized turbulence have demonstrated that submesoscales have a significant impact on mixed layer tracer budgets and hence biogeochemical dynamics, including and perhaps especially under wind forcing (e.g., Brannigan, 2016;Franks & Walstad, 1997;Lévy et al., 2009;Mahadevan et al., 2008;Whitt, Taylor, et al., 2017). However, the results presented here highlight and address an important source of uncertainty affecting these previous studies: The turbulence parameterizations have neither been designed for nor validated in an ocean surface mixed layer with significant submesoscale heterogeneity.

The results presented here are an important contribution because they are derived from coupled physical-biogeochemical ocean simulations that explicitly resolve a significant range of both submesoscales and turbulent scales and thereby greatly reduce the uncertainty emanating from turbulence parameterizations. In particular, the LES presented here are all in the high-Peclet number regime, that is, Pe ≈ 50, meaning that the resolved turbulent advective vertical tracer flux is about 50 times stronger than the subgrid-scale diffusive vertical tracer flux in the mixed layer during the storm (Figure S11). Hence, the qualitative results are expected to be insensitive to further refinement of the grid resolution or changes to the subgrid-scale model. This expectation is supported by a grid refinement study (not shown), in which the first 1.7 days of the LES,NF scenario were rerun in a smaller domain (320 m wide), and the entrainment was effectively identical on grids with horizontal resolutions ranging from 2.2 to 1.1 m, which span the 1.9 m resolution used for all the simulations reported above.

The results presented here answer the scientific question set forth in the affirmative: Submesoscale heterogeneity in a frontal zone significantly modifies the physical-biogeochemical response to the storm compared to otherwise identical scenarios without the submesoscales. Most importantly, submesoscales strongly enhance the wind-driven flux of tracers from the pycnocline to the mixed layer. For example, submesoscales enhance the maximum vertical nutrient flux in the mixed layer by a factor of 2 and the turbulent diffusivity at the mixed layer base by a factor of 10. Interestingly, the decomposition of the vertical nutrient fluxes into larger submesoscales and smaller turbulent scales shows that submesoscales contribute to enhanced vertical fluxes during the storm in two ways. First, the vertical fluxes associated with submesoscale features directly transport nutrients from the pycnocline to the mixed layer, as observed in previous studies. However, these results also show that submesoscale physical variability in a frontal zone acts as a catalyst by enhancing the small turbulent-scale nutrient fluxes near and just below the mixed layer base. On the other hand, submesoscales suppress net vertical tracer fluxes after the storm by inducing a near shutdown of turbulent scale fluxes, although submesoscale fluxes actually increase after the storm. Hence, submesoscales facilitate Journal of Geophysical Research: Oceans 10.1029/2019JC015370 a more rapid emergence of vertical phytoplankton gradients within the mixed layer after the storm. Consequently, mean phytoplankton accumulation rates are roughly twice as fast after the storm in the scenarios with submesoscales compared to their counterparts without.

Since the column model simulations do not capture this enhanced mean phytoplankton accumulation rate in the frontal zone configuration, existing physical-biogeochemical ocean models with horizontal grid spacings of 2 km or greater and the KPP and FFH parameterizations will not effectively capture the enhanced storm-driven phytoplankton growth simulated in LES. In addition, these results demonstrate that the combination of processes added to the column model with a frontal zone, which include the destratifying EBF and vertical mixing of geostrophic momentum that would also be included in a submesoscale-permitting regional ocean model, cannot explain the enhanced fluxes and biogeochemical impacts on their own. Rather, other processes or combinations of processes must be invoked to explain the results. However, the significance of these results is far from fully understood at this point, because these results are derived from only a few case studies. In addition, since the underlying processes that explain these results are not fully understood, they may impact other properties of the upper ocean beyond biogeochemistry, including heat, salt, momentum, and kinetic energy budgets.

Comparison With Other Process Simulations

Although this is the first study to present LES that resolve substantial parts of the submesoscale range together with a fully coupled biogeochemical model (with coupled nutrient/ecosystem dynamics), the results are consistent with and build on a few previous studies that have presented LES that resolve both submesoscales and turbulence and simulate passive tracers and/or phytoplankton.

For example, the result that submesoscales enhance the lateral biogeochemical variance by an order of magnitude or 2 after the storm is qualitatively consistent with the LES simulations of passive tracers in Smith et al. (2016), where it is shown that submesoscales enhance the lateral variance of tracers initialized near the base of the mixed layer by up to a factor of 20 under weak winds and weak cooling. However, the results presented here show that submesoscales significantly enhance nutrient and phytoplankton variance during the strong wind forcing for the first time, although the submesoscale tracer variance in the mixed layer is generally lower during the storm than afterward. In addition, these results show that submesoscales indirectly enhance smaller-scale turbulent tracer variance during the storm.

In addition, the result that the submesoscales are associated with enhanced wind-driven vertical mixing near the mixed layer base is qualitatively consistent with several earlier LES in short/wide frontal zones, which have shown that strong turbulent vertical mixing of tracers extends to deeper depths than predicted by KPP (Bachman et al., 2017). However, these earlier results are attributed to symmetric instability turbulence mixing over the depth of the low potential vorticity layer H q (defined to be the deepest depth where ∫ 0 H q 𝑓 ⟨ 𝜕b 𝜕z ⟩ x,𝑦 + ⟨ 𝜕u 𝜕z ⟩ x,𝑦 ⟨ 𝜕b 𝜕𝑦 ⟩ x,𝑦 dz = 0; see also Whitt et al., 2017). In the LES simulations presented here, enhanced vertical mixing extends well below H q and even penetrates into the high stratification of the seasonal pycnocline below the mixed layer (Whitt & Taylor, 2017). In addition, these results show explicitly that the enhanced vertical fluxes with submesoscales are due in large part to enhanced turbulent scale fluxes, in addition to enhanced submesoscale vertical fluxes.

Finally, the result that submesoscales reduce net vertical fluxes after the storm is also qualitatively consistent with the conclusions of Smith et al. (2016), in which it is shown that submesoscales inhibit wind/wave driven vertical mixing of passive tracers released near the mixed layer base under weak atmospheric forcing. In addition, the results are consistent with the emerging hypothesis that submesoscales suppress net vertical fluxes under weak forcing (e.g., after the storm), but submesoscales do not suppress (as in Taylor, 2016) or enhance (as in this case) net vertical fluxes as well as smaller-scale turbulent fluxes under sufficiently strong forcing.

Comparison With Observations

It is difficult to observe all the variables and spatiotemporal scales necessary to evaluate the hypotheses emerging from these simulations, but there are many observations of the bio-optical and physical response of the upper-ocean to a summer or autumn storm from satellites (e.g., Babin et al., 2004;Carranza & Gille, 2015;Fauchereau et al., 2011;Lin et al., 2003;Lin, 2012) and a rapidly increasing but small number of subsurface bio-optical observations during and after storms from profiling subsurface floats (e.g., Chacko, 2017;Girishkumar et al., 2019;Ye et al., 2013). Like all our simulations, many of these prior observations, including Journal of Geophysical Research: Oceans 10.1029/2019JC015370 the observations of Rumyantseva et al. (2015) and Painter et al. (2016), which motivated the present study, show that autumn storms increase the surface density, deepen the mixed layer, erode subsurface chlorophyll maxima and mix subsurface chlorophyll up toward the surface (thereby increasing the fraction of column-integrated chlorophyll in the mixed layer), and trigger an increase in column-integrated chlorophyll. However, the magnitude and qualitative nature of these responses are highly variable between different storms (see, e.g., Painter et al., 2016) due to the wide variety of different physical and biogeochemical circumstances in which storms occur.

Although prior observations reveal mesoscale spatial variability in the surface temperature and chlorophyll response to autumn storms (e.g., Babin et al., 2004;Girishkumar et al., 2019;Lin, 2012;Painter et al., 2016), most observations do not provide explicit information about the mesoscale lateral density variability in conjunction with explicit information about the corresponding lateral variability of submesoscales, turbulence, and biogeochemistry during and after a storm. Without such data, it is difficult to observationally evaluate the hypotheses emerging from the LES, specifically that even relatively modest open-ocean mesoscale fronts (in comparison to areas without fronts) are associated with (1) enhanced submesoscale and turbulent-scale nutrient fluxes during a storm, (2) more intermittent turbulence and suppressed net turbulent nutrient fluxes in the mixed layer after the storm, and (3) enhanced storm-driven new production.

Nevertheless, we can consider these hypotheses in light of some rare, nearly simultaneous observations of turbulence, nutrients, and submesoscales. For example, the OSMOSIS program, which motivated these process simulations, obtained observations of turbulent nutrient fluxes and submesoscales before and after an autumn storm. In addition, the OSMOSIS program provides broader context for these focused process studies, since results are available showing the full seasonal cycle of bio-optical parameters like chlorophyll fluorescence and backscatter (Bol et al., 2018;Erickson & Thompson, 2018). The observations during the storm at the OSMOSIS site reveal highly intermittent turbulent nutrient fluxes near the mixed layer base at a single location for one day after the storm. Although Rumyantseva et al. (2015) find that the intermittent mixing is plausibly attributable to intermittent alignment between the ocean currents and the wind, our simulations with a front and submesoscales (LES,F) also reveal large spatiotemporal intermittency in turbulent nutrient fluxes near the mixed layer base, both during and after the storm even though the effective wind speed after the storm is zero. Figure 7b shows how nutrient fluxes are modulated by submesoscales during the simulated storm; qualitatively similar results are found for turbulent nutrient fluxes after the storm, but not shown; see also Whitt and Taylor (2017). In addition, the hypothesis that submesoscale variability could have contributed to the observed intermittency of the turbulent nutrient flux in the observations of Rumyantseva et al. (2015) is bolstered by satellite observations of substantial submesoscale variability in the sea surface temperature before the storm (Buckingham et al., 2017). However, subsurface observations of scales smaller than 2 km are not available.

Elsewhere, a high-resolution physical/biogeochemical survey in the vicinity of the Kuroshio front by Nagai and Clayton (2017) reveals submesoscale lateral variability in nutrient concentrations, which is qualitatively similar to the simulated submesoscale nutrient variability near the pycnocline in Figure 7a. In addition, Nagai and Clayton (2017) find that these nutrient anomalies have a rectified effect on the net turbulent nutrient flux, consistent with our LES results. However, their observations were not during/after a storm.

A future interdisciplinary observational study is needed to evaluate the hypothesis that mesoscale oceanic lateral density variability induces lateral variations in the submesoscale and turbulent-scale physical and biogeochemical response to a storm using multiscale sampling of nutrients (in various forms), turbulence, and primary productivity during the life cycle of a storm at resolutions down to ∼1 km over a region ≲100 km in scale. Such a data set is challenging to obtain due to the operational challenge associated with collecting water samples during a storm, but it may be possible with an array of robotic towed and/or autonomous samplers.

Remaining Uncertainties Associated With the Model Configurations

Finally, it is worth reiterating that the results of these simulations should be viewed with caution due to a variety of outstanding uncertainties that have not been addressed here. Perhaps most significantly, the physics may be influenced by the small domain size and the absence of larger submesoscale and mesoscale structures. In addition, the biogeochemical model is at best a qualitative representation of the real biogeochemical dynamics, and it is expected that the biogeochemical model may influence the details of the Journal of Geophysical Research: Oceans 10.1029/2019JC015370 results. Further, we have only one initial condition and two fairly idealized forcing scenarios, which omit surface wave effects.

If general results are to be obtained, then future work must explore sensitivities to the details of the forcing, including wind stress direction and frequency content, as well as air-sea buoyancy flux and surface wave effects. With regard to waves, real storms are expected to generate strong wind waves, which substantially modify the Lagrangian mean shear in the upper ocean. Hence, as several previous studies have shown, waves can substantially modify the upper-ocean turbulence (e.g., Li et al., 2012;McWilliams & Sullivan, 2000), submesoscale/turbulence interactions (Hamlington et al., 2014;Haney et al., 2015;Smith et al., 2016;Sullivan et al., 2019), and thereby vertical tracer transport. It is difficult to say how waves might impact the results in this scenario, because no previous study has considered the impacts of waves on submesoscale/turbulence interactions under strong storm winds, as in this scenario. However, we expect that the addition of wind waves, roughly aligned with the winds, would be an additional source of energy for smaller-scale turbulence and enhance entrainment. But it is less clear how these waves would impact the relatively strong flow of kinetic energy to submesoscales during the storm in these simulations. Hence, a high priority for future work should be to explore submesoscale/turbulence interactions with strong waves and strong winds, building on studies with weak wind (e.g., Hamlington et al., 2014;Smith et al., 2016).

In addition, future work will have to explore the sensitivity of the results to the strength of the initial horizontal and vertical buoyancy gradients, the initial submesoscale perturbations, and the domain size and the presence/absence of larger scale structures. Future work should also explore the sensitivity to the biogeochemical parameters and model formulation.

In summary, the results presented above demonstrate that submesoscale heterogeneity can modify the mixing-induced increase in potential energy and thereby the biogeochemical response to a storm by order-one fractions. These results are novel and because they are from explicit coupled physical-biogeochemical LESs, which are the first to couple an NPZD model to LES. In addition, comparisons between the LES and a vertical column ocean model show for the first time that submesoscale/ turbulence interactions, which are not captured by current theory and ocean model parameterizations but simulated by LES, are responsible for the enhanced storm-driven biogeochemical response. Thus, submesoscale/turbulence interactions that are not represented in ocean/Earth system models may modify phytoplankton seasonal phenology in storm-driven midlatitude oceans. Since the biogeochemical impacts of submesoscale/turbulent scale physics in these idealized experiments is rather strong, future work is needed to assess the broader significance, understand the sensitivity of the results to various parameters and values not considered here, and parameterize these cross-scale physical interactions as well as the combined physical-biogeochemical interactions.

Thus, further experiments to address outstanding uncertainties with expensive submesoscale-andturbulence resolving BLES are justified. Nevertheless, a thorough exploration of these sensitivities using LES may not be feasible. A promising alternative direction is to use high-resolution LES like those presented here selectively to facilitate the development and validation of ocean models with horizontal resolutions between 10 and 100 m, which can be used to broadly explore the parameter space at a reasonable cost if they perform well in comparisons with LES. However, even this latter approach may require unconventional and expensive nonhydrostatic physical/biogeochemical ocean models. Nevertheless, nonhydrostatic ocean models at 10-100 m resolution are 2-4 orders of magnitude cheaper to run than nonhydrostatic LES at ∼1 m resolution in the same domain. In addition, some studies have used hydrostatic ocean models with horizontal resolutions down to about 100 m, although the validity of hydrostatic ocean model solutions of the ocean boundary layer at the scales below 2 km is not well understood. Thus, an important outstanding question is if and to what degree a hydrostatic ocean model with ∼100 m resolution, and therefore only a 10 by 10 or 20 by 20 grid over the 2 km box and a turbulence parameterization, can represent the dynamics simulated using BLES in this paper.

Journal of Geophysical Research: Oceans

10.1029/2019JC015370 parameter set in Table 2 was ad hoc; parameters are constrained to be reasonably consistent with published data, prior model studies focusing on the eastern North Atlantic, and via inspection of the resulting equilibrium modeled profiles in comparison with the observations reported in Rumyantseva et al. (2015) and Painter et al. (2016) and other available observations. The final parameter set in Table 2 was chosen primarily by identifying (by eye) equilibrium vertical profiles of nutrient (N) that matched the observed dissolved silica (DIS) profile reasonably well (Figure 3). The chosen parameters yield an initial profile that achieves this goal. Although more rigorous optimization is certainly possible, such an effort is beyond the scope of this paper and in our opinion is unlikely to yield qualitatively different results in the experiments presented here. The results are derived from comparisons between different physical models with the same biogeochemical model and do not depend on the model being able to exactly reproduce observations as long as it is qualitatively reasonable. A rationale for some of the chosen parameters based on available data and previously published literature is reported below.

The surface photosynthetically available radiation (PAR) is set to I 0 = 50 W/m 2 , which is equal to the climatalogical PAR from the pixel nearest the study site (48.7 • N, 16.4 • W) averaged between September and October as obtained from the VIIRS instrument (Ocean Ecology Laboratory, Ocean Biology Processing Group NASA Goddard Space Flight Center, 2018). The e-folding depth scale for the attenuation of light with depth is set to a constant 1∕k w = 18.2 m, which is intermediate between the observed e-folding depth of about 10 m (Rumyantseva et al., 2015) and the e-folding depth in clear water in the region of about 30 m (Fasham & Evans, 1995). For the parameters in Table 2, the compensation depth where phytoplankton photosynthesis is equal to phytoplankton cellular respiration, that is, where U = 𝜎 d , and the 1% light level are both between 80 and 90 m. The initial equilibrium profiles are highly sensitive to the chosen e-folding depth 1∕k w . If the other parameters are held fixed, neither the e-folding depth chosen by Fasham and Evans (1995) nor the value observed by Rumyantseva et al. (2015) produces an equilibrium profile that is a good match to observations before the storm. For simplicity, phytoplankton have no explicit influence on the PAR profile, unlike most other NPZD models. Rumyantseva et al. (2015) observe relatively small temporal changes in the e-folding depth of the PAR profile during the storm event, so the dynamic nature of the PAR profile is thought to be a minor effect.

Thee phytoplankton parameters, 𝛼, V m , and k N , and 𝜎 d were taken from the models of Fasham and Evans (1995) and are fairly typical of NPZD-class models. The mortality parameter 𝜎 d only indirectly affects the solution (via its influence on the initial condition) in these short 6-day experiments, which are about a quarter of the phytoplankton specific mortality e-folding timescale. The zooplankton parameters, R, 𝛬, 𝛾 n , 𝜁 , and ζ are less well known than the phytoplankton parameters and are tuned somewhat to produce reasonable initial conditions, but these are still fairly typical for NPZD-class models.

The detrital remineralization depth scale (sinking velocity w d by specific remineralization rate 𝛿) in the models of Fasham and Evans (1995) ranges from 100 to 400 m. Here, a much smaller value of 2.5 m is used in order to get an equilibrium in the 80 m depth domain with such simplified forcing (i.e., with a constant vertical diffusivity and no surface mixing layer). This is mostly achieved by both reducing the sinking speed from 3-6 m/day to w d = 0.25 m/day and increasing the detrital remineralization rate from 0.01-0.06 day -1 to 0.1 day -1 . These choices help in achieving realistic nutrient and phytoplankton profiles in a one-dimensional equilibrium without seasonal nutrient replenishment by deep mixing. The chosen sinking velocity is somewhat unconventional and therefore should be viewed largely as a tuning nob that helps to create an initial equilibrium profile that is similar to observations in this study. This is a reasonable approach in this case, because remineralization time scales are generally longer than the time scales of interest here (<6 days) and neither export nor remineralization are explicitly considered in this paper.

In order to crudely represent all external physical sources of nutrient that are missing in this model, such as upwelling or vertical mixing of water with higher nutrient concentrations from deeper depths, nutrient is restored to a value of N 80 = 2.5 mmol/m 3 over a time scale 1∕𝛽 = 4 days at the bottom grid cell only, that is at 80-m depth. The objective of this restoring is primarily to achieve an equilibrium vertical nutrient profile like the observed DIS profile observed during September 2012 ( Figure 3b). The choice to tune the model parameters to fit the observed DIS profile as opposed to other macronutrients probably does not change the qualitative results very much. Observed dissolved inorganic nitrogen (DIN) and phosphorous (DIP) profiles look much the same as DIS, hence initializing with those profiles would also result in an enhanced vertical nutrient flux during the storm. However, the observed stoichiometric ratios are not perfectly linear (see Journal of Geophysical Research: Oceans 10.1029/2019JC015370 Figure 13 of Painter et al., 2016). For example, the DIN:DIS ratio is about 5:1 at small concentrations in the upper 60 m but follows a 2:1 line at higher concentrations between 60 and 120 m. In addition, the DIN:DIP ratio follows an 18:1 line. The ratios of particulate organic nitrogen (PON) to biogenic silica (BSi) are about 5:1 to 10:1 in the upper 50 m (see Figure 14 of Painter et al., 2016).

The quantitative results presented in this paper are sensitive to the chosen biogeochemical parameters, which are reflected in the one-dimensional equilibrium solution used as the initial condition for the forced experiments. However, it is beyond the scope of this paper to show explicitly how either the biogeochemical initial conditions or the results of the wind-forced experiments are sensitive to changes in the biogeochemical model parameters. The qualitative and quantitative sensitivity of one-dimensional NPZD-class models to changes in the biogeochemical model parameters is discussed elsewhere (e.g., Beckmann & Hense, 2007;Doney et al., 1996;Evans & Garçon, 1997;Franks, 2002;McGillicuddy et al., 1995); recent and relevant discussions about the sensitivity of the phytoplankton and nutrient responses to wind-driven nutrient entrainment events in a similar NPZD-class model can be found in Whitt, Taylor, et al. (2017) and Whitt, Lévy, et al. (2017;n.b., section S4.1 of Whitt, Taylor, et al., 2017). The upshot is that the quantitative biogeochemical results should not be viewed as an exact reproduction of specific ocean conditions, despite some efforts documented above to choose parameters that are relevant to the observations reported by Rumyantseva et al. (2015) and Painter et al. (2016) in particular and the eastern North Atlantic more generally. However, the modeled biogeochemical responses to the different transient forcing scenarios in different physical models, all with the same biogeochemical model equations, parameters, and initial conditions, provide perspective on the physics of the different scenarios as well as the qualitative sensitivity of the biogeochemical response to physical differences alone. Hence, the discussion of the results focuses on the qualitative features of the biogeochemical solutions and comparisons between scenarios over the precise quantitative features of the results, which can nevertheless be obtained from the figures and the public data (Whitt, 2018).