## Abstract

Within the pycnocline, where diapycnal mixing is suppressed, both the vertical movement (uplift) of isopycnal surfaces and upward motion along sloping isopycnals supply nutrients to the euphotic layer, but the relative importance of each of these mechanisms is unknown. We present a method for decomposing vertical velocity *w* into two components in a Lagrangian frame: vertical velocity along sloping isopycnal surfaces and the adiabatic vertical velocity of isopycnal surfaces . We show that , where is the isopycnal slope and is the geometric aspect ratio of the flow, and that accounts for 10%–25% of the total vertical velocity *w* for isopycnal slopes representative of the midlatitude pycnocline. We perform the decomposition of *w* in a process study model of a midlatitude eddying flow field generated with a range of isopycnal slopes. A spectral decomposition of the velocity components shows that while is the largest contributor to vertical velocity, is of comparable magnitude at horizontal scales less than about 10 km, that is, at submesoscales. Increasing the horizontal grid resolution of models is known to increase vertical velocity; this increase is disproportionately due to better resolution of , as is shown here by comparing 1- and 4-km resolution model runs. Along-isopycnal vertical transport can be an important contributor to the vertical flux of tracers, including oxygen, nutrients, and chlorophyll, although we find weak covariance between vertical velocity and nutrient anomaly in our model.

## 1. Introduction

The ocean’s vertical velocity encompasses multiple processes over a range of time and space scales, from the turbulent scale to global overturning scale. Vertical velocity is difficult to measure or estimate accurately as it is from three to four orders of magnitude smaller than the mesoscale horizontal velocity and is noisy and close to the detection limit of most sensors. Though poorly understood, vertical velocity contributes to the vertical transport of tracers that are important from a biogeochemical perspective. Several properties like oxygen, dissolved and particulate forms of carbon, and nutrients like nitrate, phosphate, and silicate exhibit strong gradients in the vertical due to depth-dependent sources and sinks and weak vertical mixing. Their upward and downward transport toward or away from the sunlit surface and air–sea interface has implications for the productivity of phytoplankton, the biological pump of carbon, the exchange of atmospheric gases, and the ventilation of the oxygen-deficient ocean interior (e.g., Calil and Richards 2010; Lévy et al. 2012a; Thomsen et al. 2016; Balwada et al. 2018).

Phytoplankton primary production depends on both light, which attenuates exponentially with depth, and nutrients, which are abundant at depth and depleted in the surface sunlit (euphotic) layers. In nutrient-limited (oligotrophic) regions, phytoplankton depend on physical mechanisms to transport nutrients from depth into the sunlit euphotic zone. But, the mechanisms for vertical transport of nutrients have long been debated. Beneath the surface mixed layer, and in the pycnocline, diapycnal nutrient flux is weak, and the flow is to a good approximation adiabatic (isopycnal following) and incompressible (Lewis et al. 1986; McGillicuddy et al. 1998). In the stratified pycnocline, the vertical velocity is largely a combination of vertical motion of isopycnal surfaces (isopycnal uplift) and vertical motion along sloping isopycnal surfaces. Uncovering the mechanisms for nitrate transport in this environment could help reconcile a paradox in subtropical gyres, where carbon-based measurements of phytoplankton productivity are higher than nitrogen-based estimates, suggesting a missing nutrient source (Zhou et al. 2013).

On large scales, dissolved nitrate increases with density and depth in the ocean; the variance of nitrate is lower on isopycnal surfaces than on depth surfaces (Omand and Mahadevan 2015). Deeper isopycnal surfaces, which have higher nitrate, are also iso-nutrient surfaces. Since measured diapycnal diffusion is too weak to supply the nutrients required for phytoplankton production (Caldwell 1983; Ledwell et al. 1993), eddies and ocean fronts are thought to play an important role for the vertical nutrient flux. One mechanism for nutrient supply is the heaving of isopycnal surfaces upward into the euphotic zone, such as by mesoscale eddies. This mechanism, referred to here as “isopycnal uplift,” and previously termed “eddy pumping” (McGillicuddy et al. 1998; Charria et al. 2003; McGillicuddy 2016), is hypothesized to be of primary importance. While originally formulated at the scale of mesoscale eddies and modeled using quasigeostrophic equations (McGillicuddy et al. 1998), isopycnal uplift can also operate at smaller scales than the Rossby radius in quasi-balanced flows (Johnson et al. 2010; Ascani et al. 2013). The rate of uplift compared to the rate of nutrient uptake and the nonlinearity of the eddy determine the total nutrient transport into the euphotic zone due to isopycnal uplift (Martin and Pondaven 2003). Internal waves and tides also contribute to the vertical movement of isopycnals and can affect phytoplankton productivity (Denman and Gargett 1983; Flierl and McGillicuddy 2002), but can be separated from quasi-balanced motions associated with mesoscale eddies (Flierl and McGillicuddy 2002). They are not considered in this study, as their time scales are shorter than the typical time scales of phytoplankton growth in oligotrophic regions.

An alternate hypothesis is that transport along sloping isopycnal surfaces supplies nutrients to the euphotic zone. Though geostrophic flow on isopycnal surfaces is along pressure contours, deviations from geostrophy allow the flow to transport nutrients up- or downslope along sloping isopycnals. The degree of vertical flow is associated with the loss of geostrophic balance, and hence associated with Rossby number and submesoscale dynamics that is intensified in the vicinity of strongly sloping isopycnals or fronts and described by primitive equation models (Mahadevan and Archer 2000). Since the uptake of nutrients on an isopycnal is depth dependent, gradients in nitrate are generated on sloping isopycnals and facilitate downgradient (upslope) nutrient transport into the euphotic layer. Recent observational and modeling studies have suggested that episodic nutrient input associated with submesoscale processes could be the source of “missing” nutrients in oligotrophic regimes (Calil and Richards 2010; Ascani et al. 2013).

Much of the early work that demonstrated the importance of submesoscale physics for biogeochemical tracer flux used model simulations with progressively finer grid resolution in order to resolve smaller-scale physical processes. These studies found that productivity can be increased by resolving submesoscale physics (e.g., Mahadevan and Archer 2000; Lévy et al. 2001). Subsequent work has shown that increasing model resolution may not necessarily enhance phytoplankton productivity, since, on long time scales, increased vertical transfer of nutrients can result in a deeper nutricline (Lévy et al. 2012b) and submesoscale vertical motion may also subduct unconsumed nutrients or phytoplankton out of the euphotic zone, particularly in upwelling regions of the ocean (Gruber et al. 2011; Lathuiliere et al. 2010).

In quasigeostrophic turbulence, the vertical velocity *w* scales with the aspect ratio (of depth *H* to length *L* scale) and the Rossby number (Ro) as , and is therefore, orders of magnitude smaller than the horizontal velocity *U*. However, in dynamically unbalanced or submesoscale flows, velocities of up to 100 m day^{−1} have been observed (D’Asaro 2001; D’Asaro et al. 2018) and modeled (Mahadevan and Tandon 2006; Capet et al. 2008; Klein and Lapeyre 2009) where . While mesoscale motion is important for phytoplankton productivity, submesoscale motion could be equally important due to the much shorter time scales of vertical nutrient transport that are better matched with the phytoplankton growth rate (Mahadevan 2016).

Isopycnal uplift and along-isopycnal transport may occur in conjunction. The localized uplift of an isopycnal layer into the euphotic layer is associated with a concurrent squeezing and spreading apart of isopycnals to accommodate the motion if the motion is baroclinic. This generates a divergence/convergence of the flow along isopycnal surfaces and thereby along-isopycnal vertical motion. A water parcel uplifted into the euphotic layer on an isopycnal will be rapidly depleted of nutrients, but nutrients must be resupplied along the sloping isopycnal surface to maintain biological productivity, especially if mesoscale eddies continue to be productive for weeks. Advection along isopycnal surfaces is known to be important on the scale of an ocean basin (Palter et al. 2005; Pasquero 2005; Oschlies 2002), but it also plays a role in resupplying mesoscale features, which has not yet been quantified.

A mechanistic decomposition of the uplift and along-isopycnal vertical motions would be beneficial in improving our understanding of these processes and their representation in numerical models. It is not clear whether these kinds of vertical motion operate at different scales or in different dynamical regimes, given that isopycnal uplift can also occur on scales much smaller than the internal Rossby radius and that along-isopycnal transport can resupply nutrients to mesoscale eddies. A spectral analysis of each individual kind of motion could then reveal the scales at which each is most efficient. Furthermore, identifying their contributions to the transport of nutrients in the ocean will enable us to better quantify the physical nutrient supply. In the longer term, this will also help us design strategies to measure the vertical flux of nutrients or other biological properties in the ocean in addition to improving model parameterizations.

We begin, in section 2, by presenting a method for decomposing the vertical velocity *w* into two components, isopycnal uplift (which could act upward or downward), and vertical motion along-isopycnal surfaces . The subtlety is that these motions are so convolved that they are separable only in the Lagrangian frame, and furthermore, could occur in opposite directions with cancellation between the two. In section 2b, we derive a scaling estimate for as a fraction of the total vertical velocity *w*. In section 3, we describe a setup for the Process Study Ocean Model (PSOM) that we use for simulating the eddy field in an oligotrophic region with horizontal grid resolutions of 4 and 1 km. By applying our decomposition to the numerical model output within the pycnocline, we are able to distinguish the contribution of each vertical velocity component. We couple the physical model to a nutrient–phytoplankton model and calculate the vertical nutrient flux and resulting phytoplankton production associated with each of the vertical velocity components. Wavenumber spectra and probability density functions (PDFs) are used to characterize the vertical velocity components (section 4). Finally, in section 5, we discuss the scaling for in the global context, along with limitations and caveats, before concluding in section 6.

## 2. Vertical velocity components

### a. Decomposition

Vertical velocity can be decomposed into two components: the movement of isopycnal surfaces , which generates a vertical flux of density and of passive tracers, and vertical motion along isopycnal surfaces , which only contributes to the flux of passive tracers. This calculation must be done in a Lagrangian frame in order to separate the horizontal translation of features, such as nonlinear eddies, from Lagrangian vertical motion. In the Eulerian frame, the translation of an eddy past a point would mistakenly indicate isopycnal uplift or downdraft, though water parcels circumnavigating an eddy may not experience any vertical motion. In this study, we use the rotational (or nondivergent) horizontal velocity field as the Lagrangian frame. A velocity field which is purely rotational does not have any vertical motion and the rotational component of the flow field accounts for the bulk of the horizontal translation of features. In another setting, a different Lagrangian frame might be a more appropriate choice.

The horizontal velocity can be separated into a rotational (nondivergent) component **u**_{rot} and a divergent component **u**_{div} using a Helmholtz decomposition. To calculate the vertical motion along isopycnal surfaces, we assume that the flow is adiabatic (i.e., does not cross isopycnals). If is the height of the isopycnal *σ* surface, its horizontal gradient is and the angle the isopycnal surface makes with the horizontal *x* axis and *y* axis is given as .

The vertical component of the along-isopycnal velocity is calculated in two steps. In the first step, we project the divergent part of the total horizontal velocity along with the total vertical velocity *w* onto the isopycnal surface (Fig. 1) to obtain the along-isopycnal velocity vector **u**^{σ} = (, ) as

In the second step, the vertical component of the along-isopycnal divergent velocity is obtained by projecting the vector in the direction of the steepest slope of the isopycnal surface (which makes an angle with the horizontal) and calculating its vertical component as

The total vertical velocity is

Since ,

which simplifies to

Equations (5) and (7) enable us to diagnose and from the horizontal velocities and the slope and movement of the isopycnal surfaces. A caveat about this decomposition is that and are not mutually independent and there are situations where the two could act in opposite directions. The choice of the Lagrangian frame does not guarantee orthogonality, and admittedly, the divergent component of the velocity may also advect isopycnals.

### b. Scaling relationships

Scalings can be derived for to assess its importance in different flow regimes. To begin with, if all of the vertical motion is due to linear gravity waves, nearly all of the motion is isopycnal uplift. In past work on eddy pumping, eddies have been approximated to propagate like linear waves (Martin and Pondaven 2003). For a linear wave propagating in one dimension at constant speed *c*, the wave (isopycnal) height *h* is given by

Assuming and , we can evaluate the relative importance of uplift of the layer and vertical movement along the interface using (5) to estimate

and (7) to estimate

The ratio between the two vertical velocity components is then

Waves moving faster than individual fluid elements, , are linear. As is commonly assumed, the movement of the fluid elements can be neglected when the slope of the wave so the along-isopycnal vertical transport is negligible.

To derive a more general scaling for the relationship between the along-isopycnal component of the vertical velocity and the total vertical velocity, we assume that the flow is in quasigeostrophic balance.

First, we can estimate the magnitude of the ageostrophic velocity from the inviscid, quasigeostrophic *x*-momentum equation

where is the geostrophic velocity. The ageostrophic velocity approximates the divergent velocity to the extent that the ageostrophic velocity is irrotational. Assuming that the flow is steady and in thermal wind balance, . Here is the lateral buoyancy gradient in keeping with the notation for the vertical buoyancy gradient, , where . This implies

where is the Rossby number. The scaling for is then

where is the slope of the isopycnals. A scaling for the total vertical velocity can be derived from the continuity equation , as

which, using the scaling for the geostrophic horizontal velocity, is

The ratio is then

For a further discussion of the physical relevance of this scaling, see section 5b. In the next two sections, we evaluate this scaling relationship and investigate the role of along-isopycnal vertical motion for nutrient fluxes using a numerical model.

## 3. Methods

We use the PSOM (Mahadevan et al. 1996a,b) to simulate an eddying flow field that is analyzed for the vertical transport of nutrients into the euphotic zone. The numerical model is initialized and forced with a range of horizontal buoyancy gradients to generate fronts and eddies of varying strength. This creates 3 scenarios that we use to test the scaling in (17).

The model is set up to represent an oligotrophic (nutrient-limited) region, where the mixed layer is shallower than the euphotic layer and most of the phytoplankton production occurs beneath the mixed layer. We couple a nutrient–phytoplankton model to the physical model to investigate how nutrients are supplied to sustain new production. The model forms a deep biomass maximum that is supplied nutrients by both isopycnal uplift and along-isopycnal transport. Growth within the biological model responds to a depth-attenuating light, which generates a strong nutrient gradient in the vertical, and along sloping isopycnal surfaces.

### a. Model setup

The numerical model is set up as a periodic channel that is 256 km in the east–west (*x*) direction with closed boundaries 320 km apart in the north–south (*y*) direction (Fig. 2b). It extends to 1000-m depth and has a stretched grid in the vertical with 64 grid cells. The vertical grid spacing ranges from 3 m at the surface to 32 m at the bottom of the domain. We perform model runs with two horizontal grid resolutions, 1 and 4 km. The along-isopycnal diffusion of tracers is modeled using the Redi scheme (Redi 1982) with a diffusivity of 1 m^{2} s^{−1} in the 1-km resolution model and 20 m^{2} s^{−1} in the 4-km resolution model. The higher resolution horizontal grid is small enough to resolve some submesoscale processes. The time step for numerical integration is 432 s.

The initial salinity and temperature profiles in the center of the domain are based on an Argo profile taken on 11 December 2015 in the subtropical gyre at 28.41°N, 67.16°W. The center of the model domain is at 35°N. The northern and southern temperature profiles are modified to create a north–south surface density gradient and two initial fronts running east–west. The initial density gradient is set with a temperature difference at the surface that decays linearly with depth to a fixed depth that varies depending on the chosen value for the surface temperature gradient. Using three different values of the north–south surface temperature variation across the 320-km domain: 2°, 3°, and 4°C, we achieve a range of isopycnal slopes (Fig. 2c) that lie within observed values. All simulations initially have a 100-m mixed layer, but isopycnal slopes differ beneath the mixed layer. The parameters used are detailed in Table 1.

Vertical motion in the ocean is dominated by internal waves, which have relatively high frequencies and large, , vertical excursions at small spatial scales. To focus on quasi-balanced vertical motion of isopycnal surfaces and vertical transport along isopycnal surfaces (section 2), we attempt to avoid generating waves in the model by initializing the flow in geostrophic and thermal wind balance, and not using any wind forcing. Once the model is spun up, we maintain a statistically steady state with a net surface heat flux (comprising incoming shortwave − outgoing longwave radiation) that varies from a positive value in the south to negative value in the north of the domain, but integrates to zero over the whole domain with no seasonal cycle and no evaporation or precipitation (Fig. 2a).

### b. Decomposition of vertical velocity

The and are calculated from the model velocity field using (5) and (7). All gradients are calculated using centered differencing. Model fields are interpolated onto isopycnals and a Helmholtz decomposition of the horizontal velocity field on a density surface separates into divergent and rotational components as

where

Solving the Poisson equation [(19)] for the velocity potential *ϕ* with periodic boundary conditions on the east–west boundaries and Neumann (no flux) boundary conditions on the closed north and south walls, gives the divergent velocity field . The vertical velocity decomposition is calculated on five isopycnal surfaces from the model output at 3.6-h intervals so , where = 3.6 h. Over the entire analysis window, there are 200 time intervals.

### c. Biological model

The biological model creates a deep biomass maximum due to a balance between light and nutrient availability, as is observed in many oligotrophic regions. The model is adapted from Hodges and Rudnick (2004) as a minimal model to create a deep biomass maximum by growth of phytoplankton *P* as a function of light and nutrient uptake (modeled using a Michaelis–Menten uptake function) as

Model parameters and values are defined in Table 2. The sinking velocity decreases linearly from at the surface () to 0 at . The quadratic phytoplankton loss term models phytoplankton concentration–dependent loss due to, for example, predation. It also increases the response to episodic nutrient inputs by reducing the phytoplankton loss rate when the phytoplankton biomass is low. The phytoplankton loss is returned to nutrients.

In this model, the euphotic depth (defined as the depth where the light level is 1% of that at the surface) is at 125 m, using a light attenuation coefficient that is typical for subtropical gyres (Dickey et al. 2001). The euphotic depth is deeper than the mixed layer and the maximum depth of the density gradient is well below the euphotic layer for all simulations. The nutrients in the model are initialized using a linear nutrient–density relationship *μ*mol kg^{−1} computed from glider time series near Bermuda (R. Curry 2016, personal communication) so that nutrient concentration is initially constant on isopycnal *ρ* surfaces. Nutrient concentration is nearly zero in the mixed layer. The biological model is initialized when the physical model is approaching equilibrium. The nutrient fluxes are analyzed once both the biological model and the physical model come into equilibrium, as diagnosed by small rates of domain-integrated change of phytoplankton concentration and kinetic energy, respectively.

## 4. Results

### a. Effects of varying isopycnal slope

Varying the isopycnal slope and the model grid resolution results in distinct flow fields (Fig. 3). The dominant eddy size is smaller in the model runs with shallower isopycnal slopes than those with steeper isopycnal slopes. The high-resolution model runs have much larger magnitudes of surface relative vorticity, attaining values as high as 2–3*f*, compared to the low-resolution model runs, which have relative vorticity of less than 0.5*f*. The high resolution model runs have more finescale vorticity filaments.

All of the simulations have very similar mean nutrient distributions during the analysis period (Fig. 4b) with deep biomass maxima just above 100 m. The simulations with larger horizontal buoyancy gradients have slightly higher nutrient concentrations at a depth of 500 m due to higher productivity and more sinking and remineralization at depth. By contrast, the simulations with weaker surface lateral buoyancy gradients have lower nutrients in the near-surface. Due to enhanced along-isopycnal diffusion in the low-resolution model runs, there is a smaller gradient with depth in the nutrient distribution, with slightly lower nutrient concentration near 500-m depth and slightly higher nutrient concentration in the surface. The model generates patchy and filamentous distributions of nutrients and phytoplankton production (Fig. 5c).

The high-resolution and low-resolution models that have the same initial condition maintain similar isopycnal slopes on the large scale. When the density is averaged in 50 km × 50 km boxes over the 30-day model run, the distribution of isopycnal slope with depth in the high-resolution and low-resolution models with the same initial condition is similar (Fig. 4a). In fact, the higher resolution models have more eddy activity and so more efficiently restratify fronts through baroclinic instability, resulting in weaker isopycnal slopes at intermediate depths. However, at the local scale, the high-resolution models develop larger isopycnal slopes and density gradients (Fig. 6).

The models with steeper isopycnal slope and higher resolution have the largest vertical velocity, as evidenced by heavier tails in the PDF of vertical velocity (Fig. 7a). As the total vertical velocity increases, so does the proportion of the vertical motion that is along isopycnal surfaces since also increases with the divergent component of the horizontal velocity and the isopycnal slope in accordance with (5).

### b. Along-isopycnal vertical velocity

In this model configuration, within a subtropical setting, we can specify the length scale as the internal Rossby radius in order to write the geometric aspect ratio and the scaling (17) as

The combination of changing the large-scale isopycnal slope and the model resolution, results in a range of values of this quasigeostrophic scaling parameter. In these simulations, the along-isopycnal vertical velocity ranges from being a negligible contribution to the total vertical velocity in the low-resolution model runs, particularly with shallow isopycnal slopes, to representing about 25% of the total vertical transport in the high-resolution model run with the largest isopycnal slope (Fig. 6). In these results, the ratio of along-isopycnal, to total, vertical velocity is domain-integrated within the pycnocline and computed on five isopycnal surfaces in the range *σ* = 25.7–26.1. The linear relationship between the model-derived ratio and (Fig. 6, slope = 2.58, ) affirms the scaling relationship (21).

The spatial structure of differs from that of (Figs. 8b,c and Figs. 5a,b). The along-isopycnal vertical velocity has a shallower wavenumber spectrum and more power at small scales (Figs. 8b,c). With a weaker isopycnal slope (Δ*T* = 2°C and Δ*T* = 3°C) the relationship between the spectrum and the spectrum is similar in the low-resolution and high-resolution simulations. With the larger buoyancy gradient, has relatively more power across all wavenumbers in the high-resolution simulation than in the low-resolution simulation. With the largest buoyancy gradients (Δ*T* = 4°C, 1-km resolution), has more power at large scales than in the other simulations and the spectra of and are nearly identical.

At the surface, the vertical velocity is skewed such that the negative vertical velocities (downward) are more intense than the positive vertical velocities (not shown). Deeper in the water column, in the depth range 100–135 m, which is near the euphotic depth, the distribution of the total vertical velocity is symmetric (Fig. 7a). However, the distribution of along-isopycnal vertical velocity is skewed even near the euphotic depth, with more intense downward vertical velocities (Fig. 7b). The skewness in is more intense in the high resolution and steep isopycnal slope simulations. Relatedly, is slightly positively skewed (Fig. 7c).

### c. Nutrient fluxes

The vertical flux of the nutrient-like tracer depends not just on the magnitude of the vertical motion, but also on the covariance between the tracer anomaly and the vertical velocity . We define the nutrient anomaly as the anomaly from the mean concentration on a given isopycnal surface () in order to focus on flux along isopycnal surfaces, and because nutrient concentration is initially correlated with density. In these simulations, the magnitude of the vertical nutrient flux induced by along-isopycnal vertical transport increases in proportion to vertical along-isopycnal motion such that , where the overbar indicates a spatial average over the domain. This is due to low coherence (calculated, but not shown) between all components of the vertical velocity (*w*, , and ) and the nutrient anomaly in this model, also evident from the large depth range of the nutricline (Fig. 5c). Here, the total nutrient flux depends on both model resolution and isopycnal slope (Figs. 7d–f). Due to the low covariance between vertical velocity and nutrient anomaly, the negative skewness in the distribution of necessarily results in a negatively skewed vertical nutrient flux due to along-isopycnal vertical velocity (Fig. 7e).

## 5. Discussion

### a. Global context

We assess the extent to which the vertical flux (of a tracer) occurs along sloping isopycnal surfaces in the pycnocline for different regions of the world’s oceans by computing the scaling factor in a high resolution global MITgcm run with 1/48° resolution in the horizontal and 90 vertical levels (Su et al. 2018). The horizontal resolution of the MITgcm is approximately 2 km at midlatitudes, which is midway between the two model resolutions used in the smaller domains analyzed in this study. While the necessity of performing a Helmholtz decomposition to calculate makes computation of the ratio computationally costly on such a large domain, the scaling factor can be easily computed from the global model fields (Fig. 9). For the global model, it is inappropriate to choose the length scale , which is based on the midlatitude baroclinic instability. Though this simplifies the scaling factor to as in Fig. 6, the simplification is applicable only in the subtropical gyres, and not near the equator and poles because of the inverse dependence on the Coriolis frequency. The inverse aspect ratio of the flow is instead calculated globally using the magnitudes of the velocity gradients as .

The extent to which the vertical motion is along sloping isopycnal surfaces varies geographically and depends on the flow regime. The scaling factor is large in the low latitudes, where the aspect ratio of the velocities becomes small. In the mid latitudes, becomes relatively large, reaching values of 40%–50% in the highly energetic regions such as the Gulf Stream (Fig. 9d) and Kuroshio, as well as the California Current and Humboldt Current (Fig. 9c). In the high latitudes, the large values of are mostly due to large isopycnal slopes. There is some seasonality to the scaling factor, which is larger during the winter, particularly in high latitudes.

This suggests that vertical motion along sloping isopycnal surfaces is particularly important in low latitudes and in highly energetic regions, or fronts, with large isopycnal slope. Though the seasonality in the scaling suggests a seasonality in the importance of along-isopycnal vertical flux, its importance for production may be diminished in winter, because nutrient flux during the wintertime is less important for productivity than in the summertime. The scaling suggests that the vertical flux along isopycnal surfaces could be the dominant vertical flux pathway in some regions.

### b. Along-isopycnal transport

Using this decomposition, we find that and are anticorrelated over much, but not all, of the domain. Although different processes drive and , it should be pointed out that they are not mutually independent, and they co-occur in some dynamical situations. For example, as the stratification changes due to an uplift process, water must move along a sloping isopycnal surface. On the other hand, can occur independently of along a steady isopycnal surface, and a linear wave is an example of a process where only is important. Though and co-occur in our model eddy field, the difference in the distribution of their spatial scales suggests process-level differences in their contribution to vertical motion.

Much of the subtropical gyres are nutrient-limited and it has been proposed that the eddy-driven supply of nutrients is particularly important in these regions (McGillicuddy et al. 1998; Ascani et al. 2013). In the case of midlatitude baroclinic instability, using the thermal wind relation and , enables us to express (17) in terms of , the thermally balanced Richardson number as

With steeper isopycnal slope, the total vertical velocity increases, but so does the proportion of the vertical velocity that is due to motion along sloping isopycnal surfaces. In our model runs, the largest value of , equivalent to . This raises questions about the appropriateness of this quasigeostrophic scaling at smaller Richardson number and in higher resolution models, though it may hold beyond the quasigeostrophic regime. In the semigeostrophic solution of the ageostrophic overturning circulation, which is more appropriate for small Richardson number, the flow is more aligned with isopycnal surfaces than in the quasigeostrophic regime (Hoskins 1982). This suggests that more of the vertical flux is along-isopycnal, rather than due to the motion of the isopycnal surfaces. If we assume that the equilibrated state of mixed layer restratification by inviscid geostrophic adjustment is (Tandon and Garrett 1994), we find that , even as . However, the scaling relationship does not account for processes like mixed layer instability, which can restratify a weakly stratified mixed layer more efficiently than geostrophic adjustment (Fox-Kemper et al. 2008). Other processes, such as near-inertial oscillations (Tandon and Garrett 1995), wind-induced Ekman transport (Whitt et al. 2017), and symmetric instability (Thomas et al. 2013) may also contribute to the along-isopycnal nutrient fluxes, but are not considered here. Moreover, mixed layer processes could be important for vertical tracer flux even beneath the mixed layer through the coupling of submesoscale and mesoscale processes (Ramachandran et al. 2014).

### c. Spatial scales

Previous studies (Mahadevan and Archer 2000; Lévy et al. 2001; Rosso et al. 2014) have shown that the magnitude of the vertical velocity and vertical nutrient flux increases with increasing model resolution. Our results are consistent, but additionally reveal that is more poorly represented in low-resolution model simulations than . Since has a relatively flat wavenumber spectrum compared to that of , representation of small scales in models is particularly important for resolving the contribution to vertical velocity from . The high-resolution model runs generate larger density gradients at the local scale, resulting in a higher proportion of the vertical motion being due to . In addition, in the simulation with the largest horizontal buoyancy gradient, representation of the smallest spatial scales results in more power of across spatial scales.

### d. Nutrient fluxes

These results suggest a path forward for parameterizing subgrid scale tracer flux in mesoscale eddy resolving and eddy permitting models that is distinct from the parameterization of buoyancy fluxes. The Redi scheme (Redi 1982), which is used in our model runs, parameterizes horizontal tracer mixing as diffusion acting along isopycnal surfaces, with the vertical diffusion of tracer proportional to the isopycnal slope squared (Gnanadesikan et al. 2015). We find that in addition to the isopycnal slope, the velocity gradients are also important for determining the vertical tracer flux, although these may be less important for horizontal tracer fluxes. Parameterizing the along-isopycnal tracer flux is important for appropriately representing the isopycnal uplift flux as well, because both components of the vertical motion contribute to the nutrient anomalies.

In these model simulations, the coherence between vertical velocity and nutrient anomaly was relatively low. The nutrient anomalies may persist for longer times than the vertical velocity anomalies, which have high-frequency fluctuations (Pasquero 2005). For this reason, the coherence between vertical velocity and nutrient anomalies may depend on the nutrient uptake time scales, with faster uptake resulting in a higher coherence between the vertical velocity and nutrient anomaly. The coherence between the vertical velocity and the nutrient anomaly is very noisy in the high-resolution simulation, due to a very filamentous nutrient distribution along isopycnal surfaces that is not necessarily in phase with the vertical velocity. In the case of a tracer-like nutrient, an increased rate of uptake will enhance the vertical gradient of the tracer and lead to more coherence between the vertical velocity and nutrient anomalies.

Along-isopycnal transport of nutrients could be particularly important for supplying deep biomass maxima, which are common in oligotrophic regimes and tend to occur around 100-m depth (Cullen 2015). The depth of the mixed layer relative to the depth of light availability and the nutrient content of the mixed layer are known to be important influences on phytoplankton growth (Sverdrup 1953). Since production can be enhanced by along-isopycnal transport, the depth of lateral density gradients relative to the depth of the nutricline may be another important factor for total productivity.

In these simulations, particularly those with the largest density gradients, the nutrient flux is negatively skewed due to the negatively skewed vertical velocity field combined with a low coherence between nutrient anomaly and vertical velocity. Some previous studies have found that increasing grid resolution in highly energetic regions can decrease phytoplankton productivity (Gruber et al. 2011; Lathuiliere et al. 2010) due to an export of phytoplankton or unutilized nutrients. This is likely due to contributing to a negatively skewed flux and the along-isopycnal export of nutrients before the nutrient can be consumed.

As the measurement of vertical velocity in the ocean becomes possible through Lagrangian floats (D’Asaro et al. 2018) or direct measurement (Thurnherr et al. 2015), this framework could be applied to observational data to separate vertical transport due to linear waves and wavelike vertical motion by eddies, from tracer transport along sloping isopycnal surfaces. Observations by profiling floats have been able to observe motion of isopycnal surfaces and motion of nutrient isosurfaces (Ascani et al. 2013). While these may be correlated, the link is not necessarily causal. Instead, the vertical flux of tracers along sloping isopycnal surfaces should be understood as potentially distinct from the vertical flux of buoyancy.

## 6. Conclusions

The decomposition of the subsurface vertical velocity, *w*, into an along-isopycnal component , and an isopycnal uplift component , enables us to mechanistically understand vertical motion and derive a scaling for the ratio . Numerical simulations performed with 1- and 4-km grid resolutions satisfy the scaling relationship for a range of isopycnal slopes. The scaling suggests a new parameterization for the vertical component of tracer transport along isopycnal surfaces in eddy-permitting and eddy-resolving models. accounts for 5%–25% of the vertical transport of nutrients in simulations of the midlatitude eddy field, but is potentially the dominant mechanism for vertical transport in regions with steeply sloping isopycnals. While steeper isopycnal slopes and higher model resolution results in larger, and more negatively skewed vertical velocity, *w*, we find this is mostly due to the increased magnitude and negative skewness in compared to . Vertical motion along sloping isopycnal surfaces is particularly important at small scales, in regions with large lateral density gradients, and at low latitudes. While does not contribute a buoyancy flux, it is important for the vertical flux of biogeochemical tracers with strong vertical gradients, and for the supply of nutrients to phytoplankton in oligotrophic regions.

## Acknowledgments

MAF was supported by a National Defense Science and Engineering Graduate Fellowship and AM by NSF OCE-I434788. The authors thank Glenn Flierl and Ruth Curry for helpful conversations, and three anonymous reviewers for comments that improved the manuscript.

## REFERENCES

*Biological-Physical Interactions in the Sea*, A. R. Robinson, J. J. McCarthy, and B. J. Rothschild, Eds.,

*The Sea—Ideas and Observations on Progress in the Study of the Seas*, Vol. 12, John Wiley and Sons, 113–185.

^{234}Th-based particle export associated with anticyclonic eddies

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).