Modulation of Phytoplankton Uptake by Mesoscale and Submesoscale Eddies in the California Current System

Eddies play a crucial role in shaping ocean dynamics by affecting material transport, and generating spatio‐temporal heterogeneity. However, how eddies at different scales modulate biogeochemical transformation rates remains an open question. Applying a multi‐scale decomposition to a numerical simulation, we investigate the respective impact of mesoscale and submesoscale eddies on nutrient transport and biogeochemical cycling in the California Current System. First, the non‐linear nature of nutrient uptake by phytoplankton results in a 50% reduction in primary production in the presence of eddies. Second, eddies shape the vertical transport of nutrients with a strong compensation between mesoscale and submesoscale. Third, the eddy effect on nutrient uptake is controlled by the covariance of temperature, nutrient and phytoplankton fluctuations caused by eddies. Our results shed new light on the tight interaction between non‐linear fluid dynamics and ecosystem processes in realistic eddy regimes, which remain largely under‐resolved by global Earth system models.


of 12
Circulation at eddy scales affects biogeochemistry in multiple ways. In the simplest way, eddy-induced physical-biogeochemical interactions occur via two main processes: eddy transport and eddy reaction rates (Levy et al., 2023;Levy & Martin, 2013). These are similar in essence, but reflect different underlying mechanisms (Goodman & Robinson, 2008). Eddy transport arises from eddy-scale correlations between fluctuations in currents and tracer concentrations. This is an advective stirring process with both vertical (Benitez-Nelson et al., 2007;Falkowski et al., 1991;Oschlies & Garcon, 1998) and horizontal (Gaube et al., 2014;Gruber et al., 2011;Lathuilière et al., 2010) contributions. The effects of eddy transport depend on the circulation regime and large-scale biogeochemical gradients, and remain an active field of investigation (Lévy et al., 2018;Levy et al., 2023).
Eddy reaction rates consist of a "rectification" of large-scale, low-frequency biogeochemical transformation rates that arises from the non-linear nature of biogeochemical interactions (which include primary production, zooplankton grazing, remineralization, and air-sea fluxes) in a turbulent, heterogeneous environment. As a result, biogeochemical transformation rates estimated from a "mean field approximation," that is, from properties averaged over scales greater than those of eddies, often fail to represent the biogeochemical dynamics of a turbulent ocean (Brentnall et al., 2003;Rovinsky et al., 1997). In analogy to eddy transport fluxes, a Reynolds decomposition can be applied to biogeochemical rates to separate mean from eddy contributions. This approach relies on appropriate spatial or temporal filters to separate the effects of the mean tracer distribution from fluctuations induced by eddies (Goodman, 2011;Goodman & Robinson, 2008;Wallhead et al., 2008).
Beyond theoretical and idealized work (Brentnall et al., 2003;Goodman & Robinson, 2008;Wallhead et al., 2008), several studies characterized the role of mesoscale eddies on processes ranging from nutrient supply and primary production (Mahadevan, 2016;McGillicuddy Jr et al., 1998), biological carbon export (Harrison et al., 2018;Resplandy et al., 2019), and air-sea CO 2 fluxes Souza et al., 2021). By using an idealized, eddy-resolving simulation of the North Atlantic Ocean, Levy and Martin (2013) showed that eddy contributions accounted for between 5% and 30% of primary production and grazing rates. Eddy effects were mostly attributed to mesoscale variability (with typical length scales of between 30 and 100 km). To date, most estimates of eddy effects on biogeochemistry have focused on mesoscale activity, while submesoscales remain under-resolved and poorly characterized. Thus, it is possible that, in region with vigorous submesoscale activity-such as intense frontal regions and upwelling systems, eddy reactions may be more important than previously appreciated.
The California Current System (CCS) is ideally suited for studies of eddy-driven physical-biogeochemical interactions. In this coastal environment, wind-driven upwelling of nutrient-rich waters fuels intense biological productivity (Carr & Kearns, 2003;Messié et al., 2009) and generates a highly energetic field of mesoscale and submesoscale eddies (Capet et al., 2008;Marchesiello et al., 2003). Baroclinic instabilities of the alongshore current (Marchesiello et al., 2003) result in a cross-shore transport of nutrients and organic material followed by subduction along the CCS fronts. This so-called "eddy quenching" mechanism (Gruber et al., 2011), which is enhanced at submesoscales , reduces productivity in the coastal band, and supplies nutrients to remote open-ocean regions (Frenger et al., 2018;Lovecchio et al., 2018;Yamamoto et al., 2018). However, the contributions of mesoscale and submesoscale eddies on transport and biogeochemical transformation rates remain poorly characterized in this region.
Here, we evaluate the role of mesoscale and submesoscale eddies on nutrient transport and phytoplankton uptake rates by applying a multi-scale Reynolds decomposition to output from a submesoscale-permitting model of the CCS (Damien et al., 2023;. Our analysis provides new insights on the routes of nutrient supply and removal in the euphotic layer, and on the scale-dependent interplay between non-linear fluid and ecosystem dynamics in a highly heterogeneous environment. As such, it offer insights for the improvement of coarse global models, and the development of novel parameterizations of eddy effects on biogeochemistry.

Physical-Biogeochemical Model
We use UCLA's Regional Ocean Modeling System (ROMS, Shchepetkin & McWilliams, 2005)) coupled online to the Biogeochemical Elemental Cycling model (BEC, Moore et al., 2004)). ROMS solves the hydrostatic primitive equations for the tree-dimensional velocity, temperature, salinity and the transport of tracers in a terrain-following coordinate system. BEC represents the 3 of 12 biogeochemical cycles of major elements (C, N, P, O, Fe, Si) resulting from the interaction of three phytoplankton and one zooplankton group.
We jointly analyze output from two simulations that cover the northern and southern U.S. West Coast respectively at 1 km resolution, sufficient to allow emergence of submesoscale dynamics . Both are obtained by dynamical downscaling of a coastwide configuration at 4 km resolution Renault et al., 2021), and share the same atmospheric forcings. A detailed discussion of the simulation set up and output is presented in  and Damien et al. (2023). Because these simulations do not include tidal forcings, the highest frequencies captured only comprise submesoscale circulation and internal waves generated within the domain. Output consists of physical and biogeochemical variables, transport fluxes, and biogeochemical rates averaged online at daily frequency.
In the model, an arbitrary biogeochemical tracer X i obeys the conservation equation: The first term on the right hand side, T(X i ) = −∇·(u X i ), represents the divergence of the advective flux, with u = (u, v, w) the velocity vector. It can be further decomposed into horizontal T h (X i ) and vertical T v (X i ) components. The second term represents vertical mixing, with κ the vertical eddy diffusivity. The third term, J i , is the sum of all biogeochemical rates that affect the tracer X i , which in turn depend on J model state variables X j .
We focus on the balance of nitrate (NO − 3 ), the main limiting nutrient in the CCS . For this variable, the net biogeochemical reaction rate is: here, J Uptk is the rate of uptake by phytoplankton, J Nit production by nitrification, and J Denit consumption by denitrification. Note that here, J Uptk is a negative rate because it removes nutrient from seawater. Therefore, it is equivalent to net primary production, but with an opposite sign, and expressed in nitrogen units. In the CCS, denitrification only occurs in the deeper parts of anoxic basins and the sediment, and is minor compared to nitrification and biological uptake. Hence, when discussing water column processes, we focus primarily on nitrification and uptake. The nitrification rate, J Nit = τ nit NO − 2 , is modeled as a linear function of nitrite (NO − 2 ) concentration, with τ nit a constant timescale. Non-linearities in nitrification arise from limitation under high irradiance in the euphotic zone, and inhibition at vanishing oxygen . Biological uptake depends on nutrients following a Michaelis-Menten kinetics and Liebig's law of the minimum, phytoplankton biomass, and a temperature-and light-dependent growth rate (see Text S2 in Supporting Information S1). Thus, uptake is highly non-linear because of the presence of bilinear (X i X j ), exponential ( ) , and hyperbolic (Michaelis-Menten) terms.

Triple Decomposition of Transport and Biogeochemical Rates
The non-linear nature of advection, nitrification, and uptake in the NO − 3 conservation equation (Equation 1) is at the root of eddy rectification effects that modulate the final rate of change of this tracer.
To separate the effects of mesoscale and submesoscale eddies, we apply a triple Reynolds decomposition based on two low-pass filters, ⋅ and ⋅ , with respective space/time scales ( ) and ( , ) (Capet et al., 2008). Accordingly, a model variable X i is decomposed into mean and fluctuating mesoscale and submesoscale components as: where By definition, ′ = 0 and ̃ ′′ = 0 . Here, ( ) and ( , ) , are chosen to separate mesoscale fluctuations from a large-scale, low-frequency mean ( ) that includes the seasonal cycle. The ′′ component represents the smallest scales and highest frequencies allowed by the model, that is, mostly submesos cales. The choice of filter scales depends on the circulation regime, and may not always perfectly separate intrinsic variability from forced motions. For example, along the U.S. West Coast, wind-driven upwelling is generally 4 of 12 considered part of the mean seasonal variability. However, short-term wind events can generate high frequency variability in circulation that overlaps with mesoscale and submesoscale motions. In our choice of filters, we were especially careful to attribute the main upwelling signal to large-scale regional variability (i.e., the mean term ) and not higher frequency fluctuations. To this end, we found a reasonable combination of temporal and spatial filter scales, defined as follows: • ( ) = (15 km, 3 months), with a centered averaging scheme, • ( , ) : (5 km, 3 days), with a centered averaging scheme.
We refer the readers to Text S1 in Supporting Information S1 that further discusses these filters and their performance, using surface temperature and vertical velocities as an illustration.
By applying these filters to model variables, biogeochemical transformation rates can be separated into mean and eddy components. For a nonlinear reaction rate J i (X j ) with dependence on multiple variables X j , j = 1, …, J and the transport divergence T(X i ), the analogous Reynolds decomposition takes the form: where the various terms are now calculated as: By adopting this filtering approach, the three terms in Equation 5 can be respectively interpreted as the contribution to the total rate caused by the large-scale mean tracer distributions (J mean and T mean ); the contribution caused by heterogeneity at the scale of mesoscale eddies (J meso and T meso ); and the contribution caused by heterogeneity at submesoscales and smaller scales captured by the model (J subm and T subm ). Specifically, the biogeochemical eddy contributions only exist as a rectification of biogeochemical rates that depend in non-linear ways on model variables. These contributions would vanish in the case of perfectly linear rates (Levy & Martin, 2013).

Amplitude and Sign of the Eddy Rectification
Assuming high frequency fluctuations of small amplitudes relative to the mean, the low frequency and large scale advective transport (T) and biogeochemical rates (J) can be approximated by a Taylor series expansion (Levy & Martin, 2013): An equivalent Taylor expansion can be written for the advection term T, leading to a typical definition of eddy transport fluxes (Capet et al., 2008). Since the fluctuations have zero mean, the first-order terms disappear. Ignoring the contribution of third-order terms, the amplitude and sign of the eddy effect depend on the curvature of the functional dependencies (encapsulated by J) and eddy correlation terms between model variables

Results
Along the CCS coast, the balance of NO − 3 in the surface layer (Equation 1) reflects a near compensation of two dominant effects: uptake by phytoplankton (J Uptk ), and supply by vertical transport (T v ) (Figure 1). Nearshore, these two terms account for about 80% of the NO − 3 balance. The mean component of J Uptk increases toward the coast (Figure 1a), reflecting high nutrient concentrations following upwelling (Figure 1d). Both mesoscale and DAMIEN ET AL. Supply of NO − 3 by vertical transport (Figures 1d-1f) shows noisier patterns, reflecting the high variability and large magnitude of advective fluxes. However, notable patterns emerge. The most significant is the positive mean NO − 3 supply (T v ) along the coastal band associated with upwelling. Submesoscale T v reduces the mean NO − 3 supply by 50%-70%. In contrast, mesoscale T v is weaker, and characterized by upwelling close to the coast, and downwelling offshore, thus reinforcing the mean vertical transport.
Based on these patterns, we distinguish between the coastal region, where nutrients are upwelled into the euphotic layers, and the offshore region, where subduction by mesoscale eddies dominates (Figure 1). This separation occurs at a distance of approximately 40 km from the coast, comparable with the width of the continental shelf (Damien et al., 2023).
In the coastal region, the main balance in the NO − 3 budget is between NO − 3 supplied by vertical advection and uptake by phytoplankton (Figure 2). Eddy terms oppose the mean uptake and nitrification rates (Figure 2): eddy J Uptk is positive and eddy J Nit is negative. The magnitude of the eddy uptake (J Uptk ) is particularly large, accounting for a reduction by ∼45% of the mean rate, and dominated by submesoscales. The largest source of NO − 3 , T v , is largely offset (−64%) by submesoscale subduction (Figure 2). Horizontal advection (T h ) is negligible, because of a balance between the mean component, which supplies NO − 3 , and the eddy component, which removes it. Offshore, horizontal transport replaces vertical advection as the main source of NO − 3 (58% of the total supply). Uptake is largely reduced compared to the coastal region, but the relative magnitudes of eddy rates remain similar.
NO − 3 delivery by the mean transport accounts for 64% of the horizontal supply, while mesoscales and submesoscales account for 26% and 10% respectively. Along the vertical direction, we observe a close balance between subduction at mesoscale and supply at submesoscale (Figure 2).
Over the course of the year, the NO − 3 balance is shaped by mean upwelling (Figure 3). The mean J Uptk and its submesoscale component show a maximum during upwelling (Figures 3a and 3b). J Nit is largely dominated by  (Figures 3c and 3d). In contrast, the mesoscale component of the reaction rates show weak seasonality, but large fluctuations on weekly timescales.
Near the coast, the mean T v and the associated submesoscale subduction are shaped by seasonal upwelling (Figure 3e). While mesoscale T v fluctuations cancel out when integrated over the annual cycle (Figure 2), they drive the total transport at weekly timescales. Offshore, seasonal variability is less pronounced. The period of maximum vertical transport follows the upwelling season (Figure 3f), when subduction by mesoscale eddies balances re-supply by submesoscale eddies.
Over the year, the horizontal NO − 3 flux from the coast to the open-ocean (Figure 3i) is largely positive. This redistribution of nutrients occurs at all scales, with a major contribution from the mean circulation (56%), followed In each panel, the red line shows the total rate (calculated online), which is equal to the sum of the three components. Units are mmol N m −2 s −1 . The light shaded area shows ± a standard deviation over the region. This is not included for the transport divergence because it is an order of magnitude larger than the regional average. Note that the y-axis of the transport divergence uses different scales on the left and right panels. Panel (i) shows the time series of the integrated horizontal

Eddy Transport
In the CCS, similar to other upwelling systems, nutrient subduction by eddies, or "quenching," plays a major role in modulating primary production (Gruber et al., 2011;Nagai et al., 2015;Renault et al., 2016). Here, we show that eddy quenching reflects two contrasting regimes: subduction of freshly-upwelled nutrients by submesoscale eddies nearshore, and by mesoscale eddies further offshore (Figures 1 and 2). Mesoscale eddies thus transport nutrients from the coast to the open-ocean, but also "bury" them below the euphotic zone (Gruber et al., 2011). Offshore, we observe a near compensation between subduction at mesoscale and delivery at submesoscale (Figure 2). This balance is evident between July and October, following upwelling ( Figure 3). As recently-upwelled nutrients travel offshore and progressively sink along isopycnals, submesoscale eddies tend to resupply them back to the euphotic layer .
Globally, submesoscale eddies have been shown to enhance both nutrient delivery to the surface (Lévy et al., 2001;Mahadevan, 2016), and nutrient and organic matter subduction in regions of strong frontal activity Omand et al., 2015) and upwelling systems Stukel et al., 2017). Here, we show that both effects coexist along a gradient of surface nutrient concentrations in the CCS. Specifically, the direction of submesoscale nutrient transport depends on the balance between biological uptake and nutrient supply. Relatively long nutrient residence times in surface layers associated with large nutrient concentrations and weak vertical gradients (as observed in nutrient-rich systems) favor nutrient removal by submesoscale currents. In contrast, short surface nutrient residence times associated with low nutrient concentrations and sharp nutriclines (typical of oligotrophic systems) favor submesoscale nutrient supply.

Eddy Reactions
In the California Current, eddies reduce the mean nutrient uptake, and thus net primary production, by about 50%. Most of this compensation (35%) occurs at submesoscale. This eddy rectification is significantly larger than suggested by previous studies, which mostly focused on the mesoscale (Levy & Martin, 2013;Martin et al., 2015). Our study is the first to directly assess the magnitude of eddy reactions using a submesoscale-permitting model and a scale-dependent separation of mesoscale and submesoscale. At coarser resolution, eddy kinetic energy is likely damped (Capet et al., 2008), thus leading to an underestimate of eddy heterogeneity and its contribution to biogeochemical rates.
Mesoscales and submesoscales are highly advective regimes that increase the variability in tracer fields, in turn responsible for a rectification of the mean biogeochemical rates. Integrated over large scales and low frequencies, eddy contributions consistently reduce the mean uptake (Figures 1 and 2). The magnitude and sign of this eddy effect result from the eddy covariance of tracers, and the functional dependencies that describe biogeochemical transformations (Equation 10, see also Levy and Martin (2013)). Because biogeochemical rates depend on several tracers in complex ways (see Text S2 in Supporting Information S1), eddy reaction rates generally involve contributions from the interaction of multiple tracer pairs. Analysis of the mesoscale contributions to NO − 3 uptake (Figure 4) shows that the dominant terms arise from the saturating response of uptake at high nutrient concentrations (Figures 4b-4d), and the negative correlation between NO − 3 and phytoplankton (Figures 4i and 4j). Specifically, the negative curvature of the Michaelis-Menten saturation function implies that, in a heterogeneous environment, high-frequency events characterized by large NO − 3 concentrations are not as important in boosting uptake, relative to low-NO − 3 events that are instead more effective at reducing it. At submesoscales, considering a time scale of the order of a day for nutrient uptake (Text S2 in Supporting Information S1), the ratio of reaction rates to high-frequency transport can be assumed to be small. Further assuming that high-frequency fluctuations occur mostly along the vertical direction z, Equation 10 can be re-stated as:

′2
(11) 9 of 12 with δz′ a small vertical fluctuation. Because the mean vertical profiles of nutrients and phytoplankton show large and opposite near-surface gradients, vertical fluctuations generate a positive submesoscale J Uptk . In contrast, the smaller amplitude of J Uptk at the mesoscale likely reflects a larger influence of horizontal rather than vertical fluctuations, where negative correlations between nutrients and phytoplankton are more ambiguous (Figure 4i).
When integrated over a seasonal cycle, we obtain strikingly consistent ratios between eddy and mean NO − 3 uptake terms (∼−0.35 for submesoscale and ∼−0.10 for mesoscale). To what extent these ratios can be generalized to different regions and circulation regimes remains an open question.

Implications
We found a remarkable compensation between mean and submesoscale terms in the net balance of NO − 3 over a seasonal cycle in the California Current System (Figure 2). This suggests that, in the productive coastal region, NO − 3 supply occurs predominately at large scales and low frequencies, while removal occurs at small scales and high frequencies. This balance is reversed offshore. While mesoscale contributions tend to cancel out over the . Cross sections, as a function of the distance from the coast and depth, of (a) the annual mean mesoscale eddy uptake, (b,e,h,k,n) the second derivative terms that modulate the (c) nutrient and (f) temperature eddy variance, (i) nutrient-phytoplankton eddy covariance, (l) nutrient-temperature eddy covariance, and (o) temperature-biomass eddy covariance at mesoscale. Following a Taylor series expansion (Equation 10, also shown at the top), the (a) mesoscale eddy uptake is approximated by the sum of second-order terms (d, g, j, m, p). Units for the uptake rate are mmol N m −3 s −1 . The thick black contour shows the nutricline, defined by a NO − 3 concentration of 1 mmol N m −3 . A companion figure comparing eddy covariance at mesoscale and submesoscale is provided in Figure S3 in Supporting Information S1. DAMIEN ET AL.
10.1029/2023GL104853 10 of 12 seasonal cycle, they generate large variability, producing extremes in both nutrient transport fluxes and uptake rates (Figures 3 and 4).
The nutrient heterogeneity caused by eddies does not necessarily promote biological productivity. Indeed, it systematically reduces it when averaged over large scales and low frequencies, thus representing a different kind of productivity "quenching" associated with non-linear ecosystem dynamics. The reasons are twofold. First, phytoplankton uptake quickly saturates at high nutrient concentrations. Second, high nutrient fluctuations are often associated with low phytoplankton biomass, which limits the potential for increased productivity. We also note that changes in productivity caused by correlations involving temperature are negligible in the open ocean, but become important along the continental margin (Figures 4g, 4m, and 4p).
More generally, we find that eddy terms are far from negligible compared to mean biogeochemical rates. This result questions the ability of coarser models to adequately represent nutrient fluxes and biogeochemical transformations. For example, a non eddy-resolving global model would likely overestimate the vertical nutrient supply and biological uptake along upwelling systems. Physical parameterizations of eddy transport (Fox-Kemper et al., 2008;Gent & McWilliams, 1990) can partially alleviate this issue. However, analogous parameterizations for eddy biogeochemical rates are in early stages of development (Wallhead et al., 2013) and are not yet applied to biogeochemical models. Historically, biases in ocean circulation have been addressed by tuning biogeochemical parameters, which thus implicitly depend not only on the choice of model equations, but also on the resolution at which models are run and evaluated against observations. Our finding of a constant ratio between eddy and mean nutrient uptake rates across a range of circulations (∼−0.35 for submesoscale and ∼−0.10 for mesoscale), and our analysis of the different contributions of tracer covariance terms to eddy rates, offer new insights for the development of eddy parameterizations of biogeochemical transformations.
The tight correlations between nutrients, temperature and phytoplankton fluctuations due to high-frequency vertical motions suggest that internal tides and waves, which we did not consider in this study, could also significantly alter biogeochemical reactions. Furthermore, we focused on biological nutrient uptake as the dominant biogeochemical transformation in the highly productive CCS. However, the dynamics of pelagic ecosystems is characterized by many non-linear processes, from food web interactions, to remineralization and microbial dynamics under low oxygen conditions, which remain untouched here. In environments naturally sensitive to multiple stressors, such ocean acidification, warming, and oxygen loss, eddy rectification of ecological processes could greatly alter ecosystem dynamics and marine habitats. Analysis of these processes requires a shift in emphasis from nutrients to carbon and oxygen balances, and from biogeochemical to ecological interactions.

Data Availability Statement
The model code used to generate the simulation is openly available in  (https:// doi.org/10.5281/zenodo.3988618). The simulations are reproducible using the setup and forcing described in Damien et al. (2023).