The Contribution of Submesoscale over Mesoscale Eddy Iron Transport in the Open Southern Ocean

Biological productivity in the Southern Ocean is limited by iron availability. Previous studies of iron supply have focused on mixed‐layer entrainment and diapycnal fluxes. However, the Southern Ocean is a region highly energetic mesoscale and submesoscale turbulence. Here we investigate the role of eddies in supplying iron to the euphotic zone, using a flat‐bottom zonally re‐entrant model, configured to represent the Antarctic Circumpolar Current region, that is coupled to a biogeochemical model with a realistic seasonal cycle. Eddies are admitted or suppressed by changing the model's horizontal resolution. We utilize cross spectral analysis and the generalized Omega equation to temporally and spatially decompose the vertical transport attributable to mesoscale and submesoscale motions. Our results suggest that the mesoscale vertical fluxes provide a first‐order pathway for transporting iron across the mixing‐layer base, where diapycnal mixing is weak, and must be included in modeling the open‐Southern Ocean iron budget.


Introduction
Ocean turbulence on scales of roughly 1-200 km is characterized by vigorous eddies, fronts, filaments, and other structures, which collectively make an important contribution to material transport. In the context of climate, it is well established that the mesoscale flows, roughly on the scales of 20-200 km, are of first-order importance to the global ocean heat budget, especially for the vertical transport of heat (Wolfe et al., 2008;Griffies et al., 2015;Liang et al., 2015), and can drive global-scale variations in ocean heat content (Liang et al., 2017;Busecke and Abernathey, 2019). In the Southern Ocean, the westerlies tend to steepen the isopycnals via Ekman pumping and make the water column more baroclinically unstable. This results in generation of mesoscale eddies, which counteract this wind-driven circulation by flattening the isopycnals. This balance plays a key role in the climate sensitivity of the global overturning circulation (Farneti et al., 2010;Zika et al., 2013;Gent, 2016;Sinha and Abernathey, 2016).
As observations and simulations are able to resolve smaller and smaller scales, the role of submesoscales, roughly on the scales of 1-20 km, in influencing the large-scale ocean stratification has become a topic of great interest. Due to the geostrophic inverse energy cascade (Kraichnan, 1967;Charney, 1971), models that resolve submesoscales tend to also have more energetic mesoscales than models with coarser resolution (Capet et al., 2008a(Capet et al., , 2008b(Capet et al., , 2008c. For instance, Lévy et al. (2010) showed that resolving the submesoscales The local versus remote question has important implications for eddy parameterization in coarse-resolution models. If the submesoscale local effect proves to be significant, we would either need to resolve the submesoscale dynamics in such models or parametrize the submesoscale effect on nutrient transport in order to get the correct tracer estimates and predictions. On the other hand, if the main impact of submesoscales is to enhance mesoscale transports via the inverse energy cascade (a remote effect), we may be able to rely on energy backscatter parametrizations which replicate the inverse energy cascade to energize the mesoscale field without explicitly resolving the smallest scales (Jansen et al., 2015;Anstey & Zanna, 2017;Bolton & Zanna, 2019;Bachman, 2019). In this paper, our aim is to quantify the impact of eddy fluxes on iron transport in the context of the Southern Ocean, using idealized submesoscale-permitting simulations of varying resolution, by carefully disentangling between local versus remote effects.
The Southern Ocean is one of the high-nutrient low-chlorophyll oceans (Field et al., 1998;Nolting et al., 1998;Moore & Abbott, 2000;Arrigo et al., 2008). Artificial fertilization experiments have shown that iron is predominantly the limiting nutrient in the open ocean region (Martin et al., 1990;De Baar et al., 1995;Coale et al., 2004;Aumont & Bopp, 2006). Recognizing the importance of iron in controlling Southern Ocean biomass and the associated biological carbon pump (Lévy et al., 2013), many past studies have attempted to quantify the iron supply by focusing primarily on boundary processes (i.e. Aeolian dust, and deposition by glacial melt, bathymetry and hydrothermal vents; Boyd & Ellwood, 2010;Boyd et al., 2012;Nishioka et al., 2013;Wadley et al., 2014). There has, however, been comparatively limited investigation into how ocean dynamics, specifically eddies and fronts, transport iron from the ocean interior, where it is high in concentration, to the surface layer. Conventional focus on ocean dynamics has been on coastal processes (McGillicuddy et al., 2015;Mack et al., 2017;Jiang et al., 2019), mixed-layer entrainment of iron due to wintertime cooling (Tagliabue et al., 2014;Llort et al., 2015Llort et al., , 2019 and storms (Carranza & Gille, 2015;Nicholson et al., 2016), mesoscale isopycnal heaving (Swart et al., 2015;Song et al., 2016Song et al., , 2018, and lateral stirring (d 'Ovidio et al., 2015;Ardyna et al., 2017Ardyna et al., , 2019. To the best of our knowledge, Rosso et al. (2014), Rosso et al. (2016) are the only studies that examine the effect of submesoscale dynamics on iron supply to the surface ocean in the open-ocean region of Southern Ocean, but with the geographical coverage limited to the Kerguelen Plateau region and no consideration of the seasonal cycle.
While it would be ideal to quantify eddy iron transport using in situ observations such as the Argo floats and GEOTRACES sections (Tagliabue et al., 2014), they lack the spatial and temporal resolution to sufficiently sample mesoscale and submesoscale features (Llort et al., 2018). Satellite observations on the other hand, while having good surface spatial coverage, cannot directly reveal vertical transport. We, therefore, turn to numerical simulations and tackle the local versus remote question by temporally and spatially decomposing the eddy "iron" fluxes. We accomplish this by applying cross spectral analysis, which allows us to examine the spatial and temporal scales of eddy transport, and the generalized Omega equation (Giordani & Planton, 2000). From the spectral perspective, mesoscale fluxes have larger spatial and longer timescales than the submesoscale. The Omega equation, in contrast, provides a dynamics-based decomposition, decomposing the eddy transport into a mesoscale component in balance with the geostrophic and ageostrophic horizontal flow and a submesoscale component associated with higher Rossby numbers (McWilliams, 2016;McWilliams et al., 2019;Chereskin et al., 2019).
In this study, we adopt the flat-bottom zonally re-entrant channel framework developed by Abernathey et al. (2011) and use spatial resolution as the parameter to modulate the eddy effects; we run three cases ranging from mesoscale to submesoscale permitting resolutions: 20, 5, and 2 km resolution. Although the configuration is a rather strong idealization, it has been successfully employed to investigate tracer transport in the Southern Ocean (Abernathey et al., 2013;Abernathey & Ferreira, 2015). The channel-like idealization is deliberate, as it lets us keep essentially the same mean flow and large-scale nutrient distribution across all simulations. This mitigates the confounding impacts of resolution on the basin-scale circulation reported in the North Atlantic study of ; in their case, primary production decreased with increased resolution largely due to a change in the mean circulation, namely the Gulf Stream separation which has been shown to be sensitive to submesoscale boundary layer processes (Renault et al., 2016;Schoonover et al., 2017). In contrast, the main potential remote effect of resolution in our simulations is to energize the mesoscale. This setup therefore provides an ideal test-bed to quantify the impacts of mesoscale versus submesoscale transport on phytoplankton ecology. In this study, our focus is on the resolution-and scale-dependence of physics and the physical drivers for eddy iron transport. A deeper examination of the ecosystem dynamics and overall controls on primary production is left to a companion paper in preparation (hereon U19).
Our paper is organized as follows. We describe the experimental setup briefly in the next section, and the physical results are shown in section 3. Detailed analysis of seasonal dynamics in the context of baroclinic instability and frontogenesis is found in section 3.1, and the generalized Omega equation we use to decompose mesoscale and submesoscale motion in section 3.2. We show the biogeochemical results in section 4 with the emphasis on vertical eddy iron transport (section 4.2), which is compared and contrasted with the vertical eddy buoyancy flux. Conclusions are given in section 5.

Physical Setup
The model setup is similar to Balwada et al. (2018), but without a topographic ridge, using the hydrostatic configuration of the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al., 1997). The channel domain (L x = 1, 000km × L y = 2, 000km × H = 2, 985m) is flat bottom and zonally re-entrant on a -plane centered around 49S (f 0 = −1.1 × 10 −4 s −1 , = 1.4 × 10 −11 m −1 s −1 ) in Cartesian coordinates. The run with 20 km resolution has 40 vertical levels with 10 m near the surface and increasing in separation to 100 m near the bottom. The 5 and 2 km runs have 76 vertical levels with 1 m near the surface and increasing in separation to 100 m near the bottom. Monthly varying sea-surface temperature (SST) relaxation with a piston velocity of 1∕3 m day −1 and zonal wind stress are applied at the surface; SST increases from 0 • C to 8 • C from south to north, and the zonal cosine-squared-shaped wind stress takes its maximum amplitude, between 0.1-0.2 N m −2 , at the center of the meridional extent and is tapered to zero at the northern and southern 50 km extent of the domain. The exact formulations are given in Appendix A1. The Leith-scheme horizontal  and vertical viscosity values of A h = 2.15 m 2 s −1 and A v = 5.6614 × 10 −4 m 2 s −1 are used. We apply no-slip boundary conditions at the channel walls and bottom, with the latter having a quadratic drag with C d = 2.1 × 10 −3 . Other parameter values are given in Table A1. The 20 km run was spun up from a state of rest for 200 years, and each run with higher resolution was subsequently spun up from the climatology of the lower-resolution run until the domain-averaged temperature field reached a statistical equilibrium.
Considering that 20 km resolution is mesoscale permitting and not resolving (grey zone), one may consider turning on the Gent-McWilliams eddy-induced velocity parametrization (GM; Gent & Mcwilliams, 1990;McDougall & McIntosh, 2001) developed to extract available potential energy (APE) due to otherwise resolved mesoscale baroclinic instability. Hallberg (2013), using a two-layer Phillips-type isopycnal model, however, showed that the GM parametrization for models in the grey zone acts to suppress the resolved eddies rather than replicate their effects. The suppression of eddy effects was effective in our case as well; turning on GM resulted in steepening of the isopycnals rather than slumping them (not shown). Since the focus of this study is to quantify the relative role of (sub)mesoscale eddy iron transport, the mixed-layer instability parametrization (MLI; Boccaletti et al., 2007;Fox-Kemper et al., 2011) is also another viable candidate. Mixed-layer instability is a type of baroclinic instability that feeds off the APE within the mixed layer and due to the depth scale and reduced stratification associated with it, MLI is inherently submesoscale. We show in U19 that the MLI parametrization, intended for isopycnal restratification and not for eddy tracer transport, does not enhance vertical tracer transport compared to the nonparametrized run. In this study, we therefore do not include results with eddy parametrizations turned on. Further details of the physical boundary conditions and parameters are given in Appendix A.1 and Table A1.
The Rossby deformation radius (R d ) in all resolution runs at the center of the domain is roughly 14 km. The radius was obtained by solving the Sturm-Liouville eigenvalue problem d dz = − 2 (2.1.1) where corresponds to the Rossby deformation wavenumber and radius is the inverse wavenumber (R d = −1 ). (z) is the vertical mode associated with each wavenumber. Equation 2.1.1 is derived from the linearized quasi-geostrophic potential vorticity equation around a state of rest and prescribing a plane-wave solution (Vallis, 2017, Sections 5.8.2, 9.4.3 in their book). We used the seasonal-zonal mean of stratification as the background state in solving equation 2.1.1. R d does not change seasonally, since the interior stratification does not have a strong seasonal cycle. The spatial resolution of 2 km is roughly an order of magnitude smaller than R d , allowing for partially resolved submesoscale dynamics (Lévy et al., 2018); relative vorticity ( ) reaches up to three times the local Coriolis parameter (f ) indicated by the Rossby number (Ro = ∕f; Figure 1d). Idealized models, such as ours, serve as a valuable tool to investigate the physical drivers of seasonality in (sub)mesoscale turbulence Brannigan et al., 2015) and their interaction with biogeochemistry (Munday et al., 2014).

Biogeochemical Setup
We couple the seasonally resolving physical simulation to the ecosystem model of Gloege et al. (2017), which simplifies the Darwin biogeochemical model (Follows et al., 2007) to the complexity of the two species ecosystem described in Dutkiewicz et al. (2009); there are two phytoplankton (diatoms and small phytoplankton) and two zooplankton functional groups respectively. The diatoms have a faster maximum growth rate but favor conditions with high nutrient concentrations, while the small phytoplankton are more resilient in low nutrient and light environments but have a slower maximum growth rate ( Figure A2). The model considers the full biogeochemical cycle of phosphate, nitrate, silicic acid, carbon, oxygen and, for the interest of our study, iron, with 31 distinct prognostic tracers advected and diffused by the flow.
Domain-wide nutrient supply (PO 4 , NO 3 , Fe, and SiO 2 ) is accomplished via a sponge layer at the northern 100 km extent of the domain. These nutrients then freely evolve in the interior following the circulation and biogeochemical cycles. The relaxation profiles for PO 4 , NO 3 , and SiO 2 were taken from the World Ocean Atlas at 45S and then interpolated onto our model vertical grid. We use the monthly climatological products down to 500 m where monthly data is available and append the annual climatology below. Monthly iron profiles were taken from the Biological Southern Ocean State Estimate (BSOSE; Verdy and Mazloff, 2017), as the GEOTRACES data set (Tagliabue et al., 2012(Tagliabue et al., , 2014 did not have sufficient temporal and spatial resolution. In an effort to compensate for the lack of dust, glacial, and bathymetric sources, we chose the vertical profile of iron at 50S of BSOSE, which had higher concentrations than at 45S, but details of the relaxation profiles ultimately did not make a difference in surface concentrations as the spun up interior iron concentration was rather insensitive to the details of the relaxation profile (not shown). The zonal mean of each product for all four nutrients were taken across 50-150E in the Kerguelen Plateau region. Photosynthetically available radiation (PAR) is prescribed at the surface as a meridional linear fit to the monthly-zonal mean of SeaWiFS product between latitudes of 45-60S taking its minimum in June and maximum in December ( Figure A2a-c).
In order to isolate the effects of iron supply by ocean dynamics from depth in the open ocean region, we did not include dust deposition, glacial melt, and bathymetric sources. Our ecosystem is consequently iron limited year round ( Figure A3) and pelagic community transition does not occur; diatoms dominate the community year round, whereas in the real Southern Ocean, silicic acid limitation likely comes into play (Moore et al., 2004; J. K. Carranza & Gille, 2015). Due to this dominance by diatoms, our springtime vertically integrated phytoplankton carbon biomass reaches its apex (timing of ⟨C p ⟩ maximum where ⟨·⟩ = ∫ · dz; Behrenfeld, 2010) in early November, roughly 1-2 months earlier than estimates from biogeochemical Argo floats in the Southern Ocean (Appendix B, Figure B1; Uchida et al., 2019). Although we believe that it would be possible to further tune the biogeochemical model, parameters tuned for lower resolution runs were not directly applicable to improve the ecosystem in higher resolution runs due to changes in surface iron concentration; vertical eddy iron transport increased, and the spring bloom tended to occur earlier in the year with resolution for the same biogeochemical parameters. Since 2 km resolution coupled to a full biogeochemical model is state of the art in terms of resolution, we decreased the growth rates from those used in Dutkiewicz et al. (2009) within the acceptable parameter range of previous studies (e.g. Dutkiewicz et al., 2009;Bennington et al., 2009;Llort et al., 2015;Gloege et al., 2017) as a one-shot attempt to achieve a reasonable spring bloom. The biogeochemical parameter values were then kept identical to the 2 km run for all coarser runs with the maximum growth rate for diatoms and small phytoplankton being 0.81 (= 1∕1.24) days −1 and 0.56 (= 1∕1.8) days −1 , respectively (Table A2). Considering the agreement in magnitude of the seasonal cycle and timing of bloom onset occurring around July ( Figure B1), we argue that our model, although idealized, serves as a valuable tool in quantifying the eddy transport of iron and interaction of physics and biogeochemistry at submesoscales, which are the focus of this study.
We spun up the biogeochemistry for additional 5 years after the dynamics had spun up until the domain averaged iron concentration reached a statistical equilibrium for the 20 and 5 km run. For the 2 km run, we spun up the biogeochemistry for the latter 2.5 years of the total spin up simultaneously with the dynamics at which the iron concentration over the meridional extent we analyze (y ∈ [600, 1, 400] km) reached statistical equilibrium. Further details of the biogeochemical boundary conditions, and parameter values are given in Appendix A.2 and Table A2.

Physical Results
We start by showing a snapshot of the local Rossby number (Ro = ∕f) in the top 300 m for the 2 km run on September 15, representative of austral winter ( Figure 1a) and February 15, as representative of austral summer ( Figure 1b). The skewness and seasonal differences of the probability density function (PDF) of Ro increase with resolution (Capet et al., 2008a;Lévy et al., 2010). Eddies are poorly resolved in the 20 km run, and thus the summer and winter PDFs for this case overlie each other. Vigorous eddies develop year-round for the 2 and 5 km runs, and in both these cases, wintertime vorticity PDFs have longer tails than their summertime counterparts (Figure 1c), indicative of the finer scale features developing during winter than

10.1029/2019MS001805
in summer. Figure 1d shows the outcropping of isopycnals extending from the north, a typical feature of the open Southern Ocean.
The effect of seasonal forcing in wind stress and temperature relaxation (equations A1 and A2) can be seen in the daily spatial means of SST and mixing-layer depth (MXLD; Figure 2). The resulting SST takes its maximum in February and minimum in September (Figure 2a), consistent with the seasonal cycle in the Antarctic Circumpolar Current (ACC) region of BSOSE (not shown). The MXLD is the depth over which strong isotropic mixing would be forced by surface wind stress and diabatic forcing; here we define this highly variable depth as the zonal 99 th percentile of the daily-averaged K-profile parametrization (KPP; Large et al., 1994) boundary layer. We argue that it is the mixing layer (MXL) and not mixed layer (ML) that is relevant for tracer subduction/obduction as MXL is the layer over which tracers are actively mixed (Balwada et al., 2018). The MXLD averaged over the meridional extent of y ∈ [600, 1400] km (in order to avoid the channel wall effects) is the deepest during September and shallowest in January ( Figure 2b) and follows roughly the seasonal cycle of SST, with a 1 month phase-shift in summer. This shows the MXL variability is driven by buoyancy rather than winds, which exhibit a bi-annual structure; the winds go through two full cycles of intensification and weakening per year (Figure 3d, Appendix A.1). The MXL becomes shallower as resolution increases, which is expected as MLI is better resolved with increased resolution (explained in further detail in the section below) and effectively restratifies the MXL. We attribute the 20 km runs having shallower MXLD than the 2 and 5 km run over the summer (January, February, and March) to the difference in vertical resolution.

Mixed-Layer Instability as the Seasonal Driver of Submesoscale Turbulence
We quantify the seasonal temporal variability by examining the kinetic energy (KE) of the system Callies et al., 2015;Rocha et al., 2016;Uchida et al., 2017). We remove the zero-th zonal wavenumber component (i.e., zonal-seasonal mean) from each snapshot output every 15 days, namely u ′ = u − u where the overbar denotes the seasonal and zonal mean. The 15 day interval was chosen as the timescale at which the autocorrelation of the daily-averaged horizontal velocity anomalies (u ′ ) at the center of the domain crossed zero (not shown); we treat each anomaly field (u ′ ) as an individual realization of the turbulence process in time. We take the zonal Fourier transform of u ′ and temporally average them to construct seasonal-mean spectra (|̂| 2 u ′ where( ·) = ∫ (·)e ikx dx is the zonal Fourier transform). Since our model is a re-entrant channel, all of our wavenumber spectra were taken in the zonal direction, without any tapering applied, using the Python package xrft (https://xrft.readthedocs.io/en/latest/) and then averaged over the meridional extent of y ∈ [600, 1, 400] km and temporally over each season. As a reference to our zonal-mean view of the Southern Ocean, the climatological zonal-wavenumber spectra of the Archiving, Validation and Interpretation of Satellite Oceanographic data (AVISO) geostrophic KE is shown as well using daily velocity fields sampled every 15 days. The AVISO zonal-wavenumber spectra were taken at latitudes between 50-60S wrapping zonally around the globe and then averaged meridionally assuming a Cartesian plane between those latitudes. The zonal wavenumber KE spectra ( Figure 3) shows three things as follows: (i) the mesoscales (O(50km)) are more energetic for higher resolution runs similar to Capet et al. (2008a, Fig. 6 in their paper), (ii) in winter, KE is higher at scales smaller than 25 km, and (iii) large scales (wavenumbers corresponding to scales above O(100km)) have the same order of magnitude as AVISO observations in the SO. Plotting the KE contained at scales smaller than 25km from the 2-km run against MXLD shows that their seasonality is in phase, i.e. high KE with deep MXL and visa versa (Figure 3c).
Following Uchida et al. (2017), we quantify the mechanism for the surface KE seasonality through baroclinic instability and frontogenesis. Frontogenesis is a process in which mesoscale stirring generates submesoscale filaments by bringing buoyancy contrasts closer, setting up localized sources for instabilities (McWilliams, 2016) and can be quantified by the frontogenesis function defined as . The frontogenesis function indicates whether the flow field increases or decreases the buoyancy gradients (Hoskins, 1982;Capet et al., 2008b;Brannigan et al., 2015), and was calculated using 15 daily snapshot outputs. Buoyancy here is a function of temperature only as we use a linear equation of state with no salinity (b = g ). root-mean square of vertical velocity for the 5 and 2 km run (we do not show the 20 km run as the amplitudes were orders of magnitude smaller). The amplitude of each quantity increases with resolution, which is expected from better resolved fronts and MLI (MLI; Boccaletti et al., 2007) with higher resolution. We see a strong seasonality in w ′ b ′ , an indicator of MLI within the MXL. In a spatial mean sense, w ′ b ′ and F s ∕|∇ h b| are both positive year round with the latter being more surface intensified. This implies that frontogenesis always acts to strengthen the buoyancy fronts on which MLI feeds off by converting APE to KE. Frontogenesis acting to strengthen the fronts over the summer when MXL is shallow is consistent with the seasonal cycle of submesoscale turbulence found by Brannigan et al. (2015). We find the largest positive values in F s , however, during winter in phase with MLI, which is at odds with Branningan's finding that the seasonality in MLI and frontogenesis are out of phase. This could potentially be due to different surface forcing conditions and/or eddy-mean flow interactions; Branningan's domain size was too small to allow for any eddy-mean flow interactions. The phasing between MLI and frontogenesis becomes clearer when we take the depth average over the top 100 m of w ′ b ′ and F s ∕|∇ h b| and shows a clear correlation with small scale KE (Figure 3c,d). This is consistent with Uchida et al. (2017), who, using outputs from a ocean-atmosphere fully coupled GCM, showed that even partially resolved MLI can modulate seasonality in surface KE.
Based on the increase in surface KE with resolution (Figure 3), we hypothesize that in the submesoscale range, MLI acts as an energy source for both (i) the inverse energy cascade, which energizes the mesoscales (O(50 km); Charney, 1971;Arbic et al., 2013;Qiu et al., 2014;Barkan et al., 2015;Callies et al., 2016) and (ii) the forward cascade, represented here by the Leith numerical dissipation scheme Pearson et al., 2017). We examine this hypothesis quantitatively by taking the seasonal and meridional mean of the zonal wavenumber cross spectra of w ′ and b ′ , and the KE spectral flux

1.2)
where u is the 15 daily snapshot output of horizontal velocity,( ·) the zonal Fourier transform, (·) * the complex conjugate and [] the real part of complex numbers. The cross spectra allows us to quantify the spatial

10.1029/2019MS001805
scales at which APE conversion to KE due to baroclinic instability is active and the KE spectral flux the direction and magnitude of KE cascade; in the framework of geostrophic turbulence, we would expect KE to cascade upscale ( < 0; Charney, 1971;Arbic et al., 2013). As was noted earlier for the zonal-wavenumber KE spectra, the zonal re-entrant configuration helps us avoid introducing artificial wavenumber modes and spurious errors by tapering, which Aluie et al. (2018) showed by comparing tapered spectral fluxes to their energetically consistent coarse-graining method. In other words, our spectral flux is exact in the zonal dimension.
The APE conversion rate in the surface 200 m takes its maxima at scales in the range 30-80 km for the 2 km run, which coincides with where seasonality in surface KE is apparent (Figures 3a, 5a,b). Associated with the surface maxima of APE conversion, there is a change in sign in the spectral flux around O(30km). This is particularly clear during winter with negative values in the spectral flux reaching into higher wavenumbers than in summer and the magnitude increases with resolution (Figure 5e-h), consistent with the findings by Capet et al. (2008c), Sasaki et al. (2014). The spatial scale at which the spectral flux changes sign shows the scale separation between the inverse ( < 0) and the forward cascade ( > 0). The 5 km run has a much lower APE conversion rate and consequently a weaker KE spectral flux (Figure 5c,d,g,h). The high values of vertical eddy buoyancy fluxes (w ′ b ′ ) in the surface layer in winter, at scales around the Rossby deformation radius (R d , equation 2.1.1), are a signature of the MLI feeding off the enhanced APE associated with deeper winter MXLD. These high values coincide with an inverse of KE cascade ( < 0) showing evidence that a large part of the KE produced by the MLI is transferred towards scales larger than R d (Figure 5a,b), resulting in better agreement between AVISO and higher resolution runs in the mesoscale range ( Figure 3a). The forward cascade ( > 0) at the smallest scales is due to the Leith-scheme viscosity . The large signal of inverse KE cascade at scales above O(100km) is likely coming from the deep mesoscale baroclinic instability (Figure 5a). We argue in section 4 that accurate representation of mesoscale dynamics is crucial for modulating the eddy iron transport.

Decomposing the (Sub)Mesoscales Using the Omega Equation
With the advent of submesoscale permitting GCMs, the relative importance of submesoscale heat flux over mesoscale has been an active topic of research; Su et al. (2018) argued that on a global scale, submesoscale vertical heat flux could dominate over the mesoscale. Although the zonal-wavenumber cross spectra suggests that vertical buoyancy flux associated with submesoscale turbulence is significant in our case as well (Figure 5a,b), we can further dynamically decompose the transport into balanced and unbalanced components. Instead of an ad-hoc temporal or spatial filter (e.g., Uchida et al., 2017;Su et al., 2018), here we use the generalized Omega equation (Giordani and Planton, 2000) to diagnose the vertical flow field in balance with the horizontal velocity field (Molemaker et al., 2010).
The Omega equation is purely diagnostic in the sense that it includes no terms with a time derivative and takes the form where Q is a function of the instantaneous horizontal velocities, buoyancy and pressure for which the exact form will be given in Appendix C. We define the inverted velocity from equation 3.2.1 as the balanced motion (w b ) and residual from the total vertical velocity as the unbalanced motion (w ub = w − w b ). A snapshot example of this decomposition at the depth of z = −211 m is shown in Fig. 6. To first order, the Omega equation behaves as a low-pass filter, as we can see from Figure 6a-c that w b captures the large-scale features in balance with the geostrophic and ageostrophic horizontal velocities (Giordani and Planton, 2000) and w ub captures the fronts with superposition of waves. One can observe spatial correlation between coherent mesoscale features in w b and local Rossby number (Ro = ∕f; Figure 6d); regions with a strong dipole feature in Ro around x ∈ [600, 700] km and y ∈ [1, 000, 1, 100] km, for example, also show strong dipole features in w b , and are associated with filamentary structures in w ub . It is no surprise that the contribution of w ub is large near the surface and bottom where MLI and boundary layer processes are active, while w b captures most of the variance in the interior (Figure 6e). As a reference, we also show the root-mean square profile of vertical velocity inverted from the quasi-geostrophic (QG) Omega equation (w qg ; Hoskins et al., 1978). We see that including higher-order ageostrophic terms captures the variance of total vertical velocity to a better extent (compare w b and w qg ). By including the effect of ageostrophic horizontal velocities, as opposed to  In order to quantify the scale separation of the balanced and unbalanced motion, it is useful to take the wavenumber power spectra of each component. From Figure 5a,b, it is apparent that wintertime has higher submesoscale activity so we will focus only on wintertime for the dynamical decomposition. Figure 7a-c show that the balanced motion has higher power at larger scales and at depths than the unbalanced motion. The unbalanced motion, on the other hand, is surface intensified. The decomposition in power spectra is not exact, as there is correlation between the balanced and unbalanced motion, that is,ŵ 2 =ŵ 2 b +ŵ 2 ub + 2ŵ bŵ * ub , but the power spectra of each component is a good qualitative indicator of the scale separation.
For the cross spectra of vertical velocity and buoyancy, however, the decomposition is exact, that is,ŵb ′ * = w bb ′ * +ŵ ubb ′ * . Note that the total vertical velocity (w) used here for the decomposition is slightly different from w ′ in Figure 5a,b, in that the seasonal mean is not subtracted out as we use the total horizontal velocity (u), potential temperature, and pressure fields in inverting for w b (equation 3.2.1); Figures 5a and 7d differ by the seasonal mean component, but we find the difference to be small as seen by comparing the two figures. The negative flux seen at the large scales above O(100km), is likely a result of there being some zonal structure in the seasonal mean. Although a temporal mean is statistically equivalent to the zonal mean due to our zonally re-entrant and flat bottom configuration, the difference indicates that the seasonal mean in the one year of data has a zonal structure. Upon decomposition, the negative flux projects onto the flux associated with w ub at depths below 400 m and at scales larger than O(100km) (Figure 7f). Since w ub here is defined as the residual from the total velocity, it includes both the unbalanced and seasonal mean motion and the negative flux suggests that the latter, namely Ekman pumping is counteracting the mesoscale eddies. The positive vertical buoyancy flux associated with the unbalanced motion is more surface intensified than the balanced, consistent with our understanding of MLI (Figure 7d

Biogeochemical Results
The main goal of our study is to examine how eddy iron transport behaves in our seasonally resolving and submesoscale permitting model coupled to a full biogeochemical model. To this end, we solve the governing Fe + ∇ · F diff , (4.0.1) where the left-hand side is the familiar form of tendency and advection of passive tracers. The first term on the right-hand side ( . Fe) is the source/sink term which comes from primary production, scavenging and complexation . Complexation is a process in which dissolved iron binds with organic ligands and has been shown to be a key process in increasing the residence time of iron in the nutrient pool (Gledhill and Buck, 2012;Völker and Tagliabue, 2015;Tagliabue et al., 2016). F diff is the diffusive flux. Since our domain is zonally re-entrant, it is natural to consider zonal-mean quantities, for which we show the zonal-seasonal mean iron budget as follows: We have no contribution from zonal advection due to the zonally re-entrant configuration, namely x (uFe) = 0. In the following analysis, we only consider the meridional extent of y ∈ [600, 1, 400] km away from the north/south walls.

The Relative Contribution of Eddy Flux in the Iron Budget
In order to put vertical eddy transport into perspective of the other terms in equation (4.0.2), we calculate the zonal-seasonal mean iron budget for winter (July, August, September) and summer (January, February, March) using daily averaged outputs (Figure 8). The eddy transport terms were obtained from daily averaged outputs by applying a Reynolds decomposition, that is v ′ Fe ′ = vFe − vFe. Figure 8a,b highlights the contribution by the vertical eddy and diffusive flux. It is clear that the vertical eddy transport (w ′ Fe ′ ; red dashed) is of first-order importance for both seasons, particularly during winter when MLI is active, in the budget. The eddy transport reaches deeper into the water column to bring up iron indicated by positive values (− z (w ′ Fe ′ ) > 0) than the diffusive fluxes ( z F diff ; blue), in our case due to KPP mixing. The diffusive flux is convergent near the surface, with KPP mixing transporting iron down the vertical gradients actively generated by the biogeochemical sink at the surface and eddy iron supply from the interior of iron. Putting these terms in perspective of the remaining terms in the budget (Figure 8c,d), the net biogeochemical source/sink term (̄. Fe) is a net sink near the surface year round due to primary production (. Fe p ) overwhelming the recy- cling source (. Fe r ; green dashed). The contribution due to horizontal eddy transport (− (v ′ Fe ′ ); red dotted) and mean advection (−∇ · (vFe); red solid) is small compared to the other terms in our simulation.

(Sub)Mesoscale Eddy Iron Transport
Given its dominant role in the budget, from here on, we focus on vertical eddy iron transport. Lévy et al. (2001), in the context of an oligotrophic ecosystem in a baroclinically unstable jet, showed that nutrient supply increased with spatial resolution of their model; vertical nutrient transport increased both along the submesoscale fronts and with mesoscale vertical velocities energized via the inverse energy cascade. Following our argument for vertical buoyancy fluxes, we also decompose the iron fluxes using the Omega equation, that is,ŵFe * =ŵ bF e * +ŵ ubF e * in order to investigate the relative roles of local versus remote effects on iron transport. The scale separation is again evident in Figure 9a-c where the unbalanced transport has its maximum at smaller scales. It is interesting to note that the vertical structure of buoyancy and iron transport by the unbalanced motions are quite different. The submesoscale buoyancy flux is more surface

10.1029/2019MS001805
intensified and therefore dominantly tied to the unbalanced motion (Figure 7c), while the submesoscale iron flux is strongest between 100 m and 300 m depth. Returning to the discussion between the local and remote effect, Figure 9d shows that, below the top 100 m (where the flux is weak), the transport associated with balanced mesoscale motion (remote) is larger than the unbalanced submesoscale component (local); however, unbalanced motions still contribute about 1∕3 of the total flux at 200 m depth.
One might question whether the large amplitude in iron transport at high wavenumbers around the depths of z = −200 m in the unbalanced motion is due to submesoscale turbulence or internal waves (Figure 9c). From the spectral perspective, we would expect transport associated with internal waves to have super-inertial frequencies. By definition, the temporal and spatial scales associated with mesoscale eddies would also be larger than the submesoscale eddies. We quantify this by taking the frequency and zonal wavenumber ( -k) spectra (Figure 10) of the eddy terms-defined as the deviation in hourly snapshot outputs from the zonal mean and seasonal linear trend (w ′ = w − w, Fe ′ = Fe − Fe). The -k spectra requires model outputs saved at high frequency for which we use hourly snapshot outputs, and buoyancy snapshot fields were only saved every 15 days. We, therefore, do not show the -k cross spectra for the buoyancy flux.
The -k power spectra of vertical velocity show differences between summer and winter at scales associated with mesoscale turbulence (O(30km), O( ∕f) ∼ 0.2) with wintertime having higher power. There is also a signal of internal waves at super-inertial frequencies (O( ∕f) > 1; Fig. 10a,b). However, when examining the cross spectra of vertical velocity and iron, which show the correlation between the two at each spatial and temporal scale, the signals at high frequency and wavenumber vanish for both winter and summer (Figure 10e,f). This implies that at z = −211 m depth which is right below the MXLD (Figure 4d-f): (i) eddy iron transport associated with mesoscale turbulence associated with scales larger than R d dominates over submesoscale, and ii) waves associated with super-inertial frequencies contribute to no net iron transport, consistent with the results by Balwada et al. (2018, Fig. 4 in their paper) which found that internal waves are inefficient in transporting passive tracers. The results were qualitatively similar at z = −180 m within the wintertime MXL where although there was an increase at larger wavenumbers, no increase at super-inertial frequencies in the vertical iron transport (not shown). Comparing the power and cross spectra of the 2 to 5 km run (Figure 10c,d,g,h), it is obvious that both the mesoscale eddies and vertical iron transport is weaker in the latter. We attribute this to the insufficient inverse energy cascade and resulting energetically weak mesoscale field as was discussed in section 3.1. Although the 5 km run shows no signal of internal waves, the signal at the Nyquist wavenumber is puzzling and could be due to grid-scale noise (Figure 10c,d) but again is not associated with any iron transport (Figure 10g,h).

Discussion and Conclusions
By running an idealized model of the Antarctic Circumpolar Current region at submesoscale permitting resolution with a seasonal cycle, our model partially resolves MLI (MLI; Boccaletti et al., 2007) and generates well formed mesoscale eddies and fronts (Figure 1). The agreement between surface KE in our model at scales larger than O(100km) with KE estimates from satellite altimetry observations improves with higher spatial resolution, as seen in the wavenumber power spectra (Figure 3a), likely as a result of a better-resolved inverse energy cascade ( Figure 5; Capet et al., 2008aCapet et al., , 2008bCapet et al., , 2008cLévy et al., 2010;Arbic et al., 2013). The domain size of our model has allowed us to partially resolve submesoscale turbulence and the inverse energy cascade associated with it. By coupling the channel model to a full biogeochemical model, we have examined the relative importance of vertical eddy iron transport associated with mesoscale and submesoscale turbulence in the open Southern Ocean, where understanding has primarily relied on vertical diffusion and MLE framework (Bowie et al., 2009;Tagliabue et al., 2012;Tagliabue et al., 2014;Llort et al., 2015Llort et al., , 2019. In order to quantify the temporal and spatial scales at which eddy transport was dominant, we took the frequency-wavenumber cross spectra of w and Fe. The spectra at depths below the wintertime maximum of MXLD ( Figure 10) showed two things as follows: (i) most of the vertical transport is at scales larger than the Rossby deformation radius and (ii) internal waves do not contribute to any transport of iron. The first point already implies that mesoscale turbulence is the dominant contributor to vertical iron transport below the MXL (Figure 8). Considering that different dynamics can have similar spatial and temporal scales, we further dynamically decomposed the eddy transport into its balanced and unbalanced component using the generalized Omega equation (Giordani & Planton, 2000) based on the assumption that mesoscale turbulence is associated with Rossby numbers smaller than unity (McWilliams, 2016;Lévy et al., 2018). At depths below the mixed layer, we found that the balanced motion accounts for more than half of the total vertical iron transport, whereas within the mixed layer where stratification is low, and MLI is active (Boccaletti et al., 2007), the relative contribution by unbalanced motion increases significantly (Figure 9d).
Although our wintertime biogeochemical consumption of iron is within the bounds of observations, it is too low during summer; Ellwood et al. (2008), Bowie et al. (2009) estimate it to be on the order of 100 mol m 2 yr −1 while it is roughly 35 mol m 2 yr −1 integrating over the top 100 m in our model (Figure 8b; green dotted line). Due to the lack of pelagic community transition, our bloom is too sharp and is insufficiently sustained over the summer (Appendix B, Figure B1). As the ecosystem is iron-limited year round (Moore, Mills, et al., 2013), one approach for increasing summer productivity in future work may be to reduce iron requirements for small phytoplankton, for example using separate iron-to-phosphate ratio per species, or allowing for a dynamical stoichiometry (Tréguer et al., 2018). Although the model of Dutkiewicz et al. (2009), configured to represent the global ecosystem, did not allow for this, it is likely that phytoplankton adaptation has occurred in the Southern Ocean where phytoplankton should be resilient to low iron concentrations (Tagliabue et al., 2014). Our sharp bloom should not affect our results qualitatively, however, as an increase in summertime biological consumption would result in larger vertical gradients of iron; assuming that mesoscale stirring is related to the background gradient of iron (w ′ Fe ′ ∼ dFe dz where here is the eddy diffusivity), mesoscale eddy iron transport would only increase.
Our findings, which emphasize the importance of eddy iron transport, are complimentary to Freilich and Mahadevan (2019, Figures 6, 7e in their paper) in which they show that isopycnal (sub)mesoscale stirring of nutrients increases with resolution. Rosso et al. (2014), Rosso et al. (2016), using an idealized biogeochemical model based on an exponential decay rate of iron, also argued for the importance of submesoscale iron transport in the Kerguelen Plateau region. Enhanced submesoscale turbulence in their model occurred due to flow-bathymetric interaction, resulting in a hot-spot of vertical iron transport downstream of the Kerguelen Islands. Away from bathymetric features, however, they showed that submesoscale eddy transport of iron was weak. Balwada et al. (2018) using a submesoscale permitting zonal re-entrant model with a topographic ridge showed that downstream of the ridge, vertical tracer transport was enhanced at higher frequency and wavenumber. Both studies imply that surface submesoscale turbulence is enhanced due to flow-bathymetric interaction. Since our model has flat bottom, we do not have geographical hot spots of submesoscale turbulence. We argue that in the zonal-mean sense away from bathymetric features in the open Southern Ocean, the eddy iron transport from the interior to the surface is dominated by both a local effect near the surface where MLI is active, and a remote effect, that is, mesoscale eddies energized through inverse energy cascade, at depths below the MXL (Figures 5, 9). In other words, our results suggest that in order to achieve accurate representation of iron fluxes, it is crucial to get the energetics of the mesoscale field right, either through resolving or parametrizing the inverse energy cascade due to mesoscale and submesoscale baroclinic instabilities. Current eddy parametrizations do not incorporate the effects of the inverse energy cascade; however, parameterizing this energy transfer is an active area of research (e.g., Kitsios et al., 2013;Jansen et al., 2015;Anstey & Zanna, 2017;Bolton & Zanna, 2019;Bachman, 2019). Due to computational contraints, GCMs at mesoscale-permitting resolutions will continue to be invaluable tools to investigate the coupling between physics and biogeochemistry on a global scale. Emphasis should be given to testing the impact of eddy parametrizations that incorporate the inverse energy cascade on tracer transport.
Based on modeling studies with higher spatial resolution (Molemaker et al., 2010;Smith et al., 2016;Brannigan et al., 2017;Balwada et al., 2018), it is obvious that even at 2 km resolution, which is state of the art when coupled to a full biogeochemical model, our turbulent field has not numerically converged; we would expect both the mesoscale and submesoscale tracer transport to further increase with resolution. It would be interesting to see whether the energetics and tracer transport change substantially as model resolution is

A1. Physical Forcing
The monthly structure of wind stress profile ( Figure A1) inspired by Sinha and Abernathey (2016) takes the mathematical form, where L y is the meridional extent of the domain, t ∈ [1, 2, ..., 12] for each month, and a = 100 km. The westerly jet takes its maximum in the meridional center of the domain and tails off to zero towards the boundaries. The sea-surface temperature (SST) is relaxed to the profile ( Figure A1), * (t, ) = where ( 0 , 1 , 2 ) = (3.25, 4.75, 8) • C and = 1.75 • C. Other physical parameters are listed in Table A1. The SST relaxation piston velocity is defined as the ratio between the relaxation timescale and vertical width of the top grid. In order to keep the piston velocity identical among runs, for the 20 km resolution runs, the relaxation timescale is 30 days with the vertical width of the top grid being 10 m and is 3 days with the grid width being 1 m for the 5 and 2 km runs.

A2. Biogeochemical Forcing
The light, temperature and nutrient limitation to the phytoplankton growth rate are implemented as  where max i , I i , T i , and N i are the maximum growth rate of phytoplankton i (diatom: i = 1, small phytoplankton: i = 2), and limitation factors by temperature, light and nutrients, respectively. The full equation for the time evolution of phytoplankton is given in Dutkiewicz et al. (2009, eqns. A1-A5 in their paper) and parameter values in Table A2. Each limitation factor takes values between unity and zero ( ∈ [0, 1]) with one meaning optimum conditions for growth. Light limitation is calculated as . Figure A2d,e shows behavior of I and max I in the range of PAR used in our model; the light limitation factors ( I ) are in the range of values presented by Dutkiewicz et al., 2009Dutkiewicz et al., . (2009 Fig. 1e in their paper. Note the difference in the units; 10 Ein m −2 d −1 corresponds roughly to 116 Ein m −2 s −1 .).
The temperature limitation is kept the same for both phytoplankton species and takes a similar formulation to the Arrhenius equation (e.g. Geider, 1987;Brown et al., 2004;Kremer et al., 2017) where T is the local temperature in Celsius and the reference temperature T ref defines the temperature at which T = T 0 . Lastly, the nutrient limitation factor is defined by the most limiting nutrient Figure A2. Surface PAR in our runs plotted against month and meridional extent (b). (c) The meridional mean of (a) and (b) with the colors of (a) and (c) correspond with each other. (d) The light limitation factor for each plankton (equations A4; Diatoms (P 1 ) in green and small phytoplankton (P 2 ) in yellow) and it multiplied by the maximum growth rate in (e). The blue line shows the difference between the two ( 1 − 2 ).
where N lim i = N N + N i with N j and N i being the concentration of nutrient j and half-saturation constant of nutrient j for phytoplankton i, respectively. Figure A3 shows that N lim iFe is always smallest among the nutrient limitation factors per phytoplankton, namely N i = N lim iFe = Fe Fe+ iFe , with a large dip during October to December when the climax of spring bloom takes place. In other words, iron is the limiting nutrient year-round, and the N factors for other nutrients are always close to unity, with only iron concentration showing a signal of the spring bloom. We also see the effect of diatoms having larger half saturation constants; N lim is lower for diatoms than small phytoplankton for each nutrient year round. Further details of the variables and notations will be left to Dutkiewicz et al. (2009), which the reader should refer to. Figure A3. Nutrient limitation factor (equation A6) for diatoms (i = 1; a) and small phytoplankton (i = 2; b) using monthly-mean nutrient outputs from the 2 km run averaged over the top 100 m. the meridional extent of y ∈ [600, 1400] km in our study, it is also interesting to note the indention in the southward progression of apex around y ∼ 300 km ( Figure B1a), which is likely due to our squared form of wind stress (equation A1); the MXL shoals towards the north and south boundaries ( Figure B1b). The MXL at the very south is deep due to year-round convection at the boundary wall. Figure B1. Hovmöller diagram of ⟨C p ⟩ for the 2 km run a, and from BGC Argo floats c. The hatch indicates grids where there were less than 10 Argo profiles over the five years of deployment by Southern Ocean Carbon and Climate Observations and Modeling and Southern Ocean and Climate Field Studies with Innovative Tools. b,d Zonal mean of surface KE and MXLD from the 2 km run and AVISO plotted against meridional axes. e Spatial and daily median of ⟨C p ⟩ for the 2 km run between y ∈ [600, 1400] km (black), 10 daily median biogeochemical Argo (green), and eight-day averaged C p outputs from CbPM (brown) between 45S-60S. Note the difference in the units. As CbPM relies on satellite Chlorophyll observations, the month of June lacks data due to poor light conditions. The green and brown shading indicated the zonal and daily interquartile range. MXLD, mixed-layer depth. OCE-1740648. We thank three anonymous reviewers who have greatly enhanced the readability of our manuscript. The model configuration is available on Github (doi:10.5281/zenodo.3266400) and simulation outputs for 15 daily snapshot and monthly averaged outputs of physical variables (v, , )