Global Warming Increases Interannual and Multidecadal Variability of Subarctic Atlantic Nutrients and Biological Production in the CESM1‐LE

Earth system models indicate that the Subarctic Atlantic Ocean ecosystems are uniquely sensitive to global warming. Nutrients, phytoplankton biomass, and phytoplankton production all decline precipitously in global warming scenarios. Superimposed on this forced response is changing internal variability that is not fully understood. Here, a large ensemble of simulations with the Community Earth System Model is used to quantify how global warming affects the interannual and multidecadal variability of phytoplankton production integrated over the Subarctic Atlantic. Surprisingly, it is found that this variability increases non‐monotonically with global warming. The increased variability of production is caused by an elevated volatility of wintertime surface nutrient concentrations, which is a consequence of a rising sensitivity of these nutrients to winter mixing and overturning fluctuations, which overcomes a reduced amplitude of these physical fluctuations with warming. Future work is needed to fully understand how internal climate variability impacts ocean ecosystems in a warming climate.

interannually and decadally in the region based on in-situ and satellite observations as well as in hindcast model simulations.Moreover, these variations are observed to be correlated with ocean-atmosphere drivers, for example, as measured by the North Atlantic Oscillation Index or the Subpolar Gyre Index (Barton et al., 2015;Drinkwater et al., 2003;Dutkiewicz et al., 2001;Hátún et al., 2009Hátún et al., , 2017;;Henson et al., 2009;Lampitt et al., 2010;Macovei et al., 2019;McKinley et al., 2018;Patara et al., 2011;Steinhoff et al., 2010;Tesdal et al., 2022;Ueyama & Monger, 2005;Williams et al., 2000).But, as stated by McKinley et al. (2018), "clear mechanistic links between biomass and standard physical variables or climate indices remain largely elusive." Due to the societal significance and high sensitivity of the Subarctic Atlantic Ocean biogeochemistry and ecosystems to global warming, it is a high priority to quantify and understand the spatiotemporal variability of Subarctic Atlantic production and nutrients and, moreover, how global warming has and will impact this variability.This information might help society to anticipate future changes and make any necessary modifications to marine resource management or ocean geoengineering plans and to better monitor and understand the causes of current and historical temporal variability, including the trends.In a step toward these societal objectives, this study quantifies how interannual (1-year) and multidecadal (20-year) variability in regionally integrated Subarctic Atlantic nutrients and biological production change with global warming in the high-emissions twenty-first-century Representation Concentration Pathway 8.5 (RCP8.5)scenario (Meinshausen et al., 2011;Riahi et al., 2011;Van Vuuren et al., 2011) executed in the Community Earth System Model version 1 (CESM1) (Hurrell et al., 2013).
Surprisingly, the integrated nutrients and production exhibit increasing interannual variability around their downward trends during global warming.This paper documents and explains these unintuitive results.

Approach
Due to the short and sparse observational record, it is challenging to quantify the interannual and multidecadal variability of ocean biogeochemistry and separate its causes, such as the forced responses to human activity and internal Earth system variations.Moreover, the characteristics of the internal variability may be changing as the global climate warms.Circumventing these difficulties, this study uses output from a large ensemble of simulations with the Community Earth System Model (CESM1-LE) (Kay et al., 2015) to quantify how the characteristics of the internal variability of ocean biological production change during the high-emissions RCP8.5 scenario.As demonstrated by several recent studies of ocean biogeochemistry (Elsworth et al., 2022;Li & Ilyina, 2018;Long et al., 2016;Lovenduski et al., 2016;McKinley et al., 2016;Schlunegger et al., 2019), a single-model large ensemble is a valuable complement to an ensemble of simulations with different models (as in the Coupled Model Intercomparison Project archive), because the differences between the simulations are attributable to internal Earth system variations without convoluting differences due to model formulation (Deser et al., 2020).

Model Experiments
The results are derived from 32 fully-coupled global Earth system simulations with the CESM at a nominal 1° resolution.The global ocean circulation is simulated by the Parallel Ocean Program (Smith et al., 2010) and ocean biogeochemistry is simulated by the Biogeochemical Elemental Cycling model (Behrenfeld et al., 2013;Long et al., 2016;Moore et al., 2013).The 32 members diverge from a single starting point due to the introduction of random small atmospheric perturbations in 1920, but this paper only uses the output of the RCP8.5 scenario .

Statistical Measures of Variability
The paper focuses on the dynamics of the interannual variability of the regionally integrated Subarctic Atlantic Ocean, which is defined by averaging over the ocean area from 45°N to 75°N and from 70°W to 10°E and averaging over each year or a subset of months within each year and binning annually.Some spatially-resolved results are provided in Supporting Information S1.
To define and separate the anomalies associated with internal variability in each ensemble member, a quadratic polynomial fit to the ensemble mean time series is defined as the forced response and subtracted from each member (Figure S1 in Supporting Information S1).To distinguish how the internal variability changes through time, the standard deviation (stdev) of the anomalies and the corresponding coefficient of variation (cov), that is, the stdev divided by the forced response, are calculated in overlapping 20-year windows centered every decade from 2020 to 2090.To distinguish interannual and multidecadal variability, the statistics are calculated on either all 640 1-year means (20 years from 32 ensemble members) or 32 20-year means.Figures S2 and S3 in Supporting Information S1 illustrate some 1-and 20-year anomaly time series and several measures of their changing spread with global warming, including the 1-and 20-year stdevs.Confidence intervals are defined by bootstrapping the ensemble members, since the data may be correlated within a run, but the members are independent.A statistically significant change in variability is identified when the bootstrap confidence interval on the ratio of the stdev centered at 2090 to 2020 excludes 1.The α values indicate the confidence level at which 1 is excluded from the interval; in percentage terms, the confidence is 100(1 − α).

Biogeochemical Variables
Biological production is quantified using the annual new production (Dugdale & Goering, 1967), that is the annually integrated nitrate tendency associated with the ecosystem integrated over the top 50 m (which is about equal to the euphotic depth; see e.g., Lee et al. (2007)).The new production closely approximates the export of particulate organic nitrogen from the top 50 m via sinking to depths below 50 m and is balanced by the physical supply of nitrate to the top 50 m from below (Eppley & Peterson, 1979).Below 50 m, the ecological nitrate tendency is generally positive, because remineralization exceeds uptake.The production is compared with the surface chlorophyll concentration averaged over the May-August growing season as a measure of phytoplankton biomass.Similar results are obtained for the net primary production and phytoplankton carbon concentration, which are not discussed.Finally, comparisons are made to the nitrate concentration averaged or integrated over the top 50 m in March, when surface nitrate is generally maximum, since the climatological new production is thought to be approximately equal to the net reduction of nitrate in the top 50 m during the growing season (Koeve, 2001;Whitt & Jansen, 2020;Williams et al., 2000).
Numerous prior studies have compared the CESM-BEC model solutions to relevant observations.For example, Moore et al. (2013) show that the mean chlorophyll concentration is distributed similarly to satellite observations, although the simulated mean May-August surface chlorophyll in the Subarctic Atlantic is about 30%-50% higher than observed (Figure S4 in Supporting Information S1).Behrenfeld et al. (2013) show that the seasonal cycle of surface chlorophyll and phytoplankton carbon is dynamically similar in satellite observations and the model.In addition, the interannual stdev of the May-August surface chlorophyll concentration is similar to or slightly greater than observed (Figure S4 in Supporting Information S1) (Fanton d'Andon et al., 2009;Maritorena & Siegel, 2005;Maritorena et al., 2010;Sathyendranath et al., 2019).In addition, Whitt and Jansen (2020) show that the mean surface nitrate concentration in the Subarctic Atlantic is realistic at the beginning of the RCP8.5 scenario.However, the vertical gradient of nitrate below the surface is somewhat higher than observed such that the nitrate concentration at 100 m is about 14 mmol/m 3 in the CESM and only about 11 mmol/m 3 in the observations (see also Moore et al., 2013).

Physical Drivers
The biogeochemical variability is interpreted by comparison with measures of the physical variability, including the March mixed layer depth (MLD), which is defined using a buoyancy gradient criterion (Large et al., 1997;Whitt et al., 2019), and the maximum of the zonally-integrated AMOC streamfunction in density space at 45°N, as in Whitt and Jansen (2020).The MLD is from March because the MLDs are deepest in March, and the deepest MLDs define the annual nutrient entrainment in the highly seasonal Subarctic Atlantic, where the regionally averaged MLD oscillates from about 20 m in summer to 300 m in winter.The annual average AMOC transport in the northward/upper limb is used because of its influence on the nutrient transport to the Subarctic seasonal thermocline from lower latitudes.
Many prior studies have compared the simulated MLD and AMOC to observations.For example, some comparisons between the simulated and observed maximum annual MLD are presented in Moore et al. (2013) and Whitt and Jansen (2020).Although there are some differences in spatial structure, the area-averaged maximum MLD is similar in observations and the model (about 300 m).An evaluation of the historically simulated overturning circulation in the CESM-LE is presented by Kim et al. (2018) (see also Whitt, 2019;Whitt and Jansen, 2020).For reference, Jackson et al. (2019) present measures of North Atlantic overturning variability from various ocean reanalyzes and Lumpkin and Speer (2007) presents an inverse estimate based on observations, both of which are useful comparative benchmarks (see also Buckley and Marshall, 2016;Cessi, 2019).The upshot is that the mean strength of overturning circulation in the CESM-LE, which exceeds 20 Sv at the beginning of the RCP8.5 scenario, is likely 10%-30% stronger than in the real ocean.On the other hand, the interannual variance in overturning is likely similar to observed, while multidecadal variance is likely somewhat weaker compared to the real ocean.

Increasing Variability of Nutrients and Production
Consistent with prior studies, the CESM-LE simulations confirm that rising greenhouse gases cause declines in the ensemble mean Subarctic Atlantic phytoplankton biomass and production as well as nutrients.For example, Figure 1 shows that between the 2020s and 2090s the growing-season surface chlorophyll declines from about 1.3 to 0.5 mg/m 3 (Figure 1a), the new production declines from about 0.8 to 0.5 mol N/m 2 /year (Figure 1d), and the March surface nitrate declines from about 13 to 8 mmol N/m 3 (Figure 1g).As described in Whitt and Jansen (2020), the declines in production and nutrients are primarily caused by the reduction in the rate of nutrient supply to the seasonal thermocline from the lower latitudes as the AMOC slows (Figure 1m).A secondary cause is the reduction in the March MLD, from about 280 m in 2020 to 160 m in 2090 (Figure 1j), which increases the fraction of the production that is lost from the seasonal thermocline via export (Palevsky & Doney, 2018).
One might expect that these strong downward trends in the biomass, production, nutrients and their physical drivers might be accompanied by reductions in the corresponding interannual and multidecadal variability (Elsworth et al., 2022).And, indeed, the interannual stdev of the chlorophyll declines by about 30% during the RCP8.5 scenario (Figure 1b) in conjunction with declines in the interannual stdevs of the March MLD and the AMOC (Figures 1k and 1n).However, the covs of the chlorophyll increase from about 6% to 11% and from about 2% to 4% at interannual and multidecadal timescales, respectively (Figure 1c).Moreover, the interannual and multidecadal stdevs of production variability actually increase non-monotonically between 2020 and 2090 from 0.029 to 0.032 mol N/m 2 and from 0.011 to 0.016 mol N/m 2 at interannual and multidecadal timescales respectively (Figure 1e).The corresponding covs roughly double from about 3.5% to 6% and 1.5% to 3% (Figure 1f).The stdevs of March nitrate increase similarly but more prominently than production between 2020 and 2090, from about 0.38 to 0.54 mmol/m 3 and 0.13 to 0.28 mmol/m 3 at interannual and multidecadal timescales respectively (Figure 1h).The corresponding covs rise from about 3% to 6.5% and from 1% to 3.5%, similar to those of production (Figure 1i).For both production and nitrate, the rate of change of the variability accelerates during the warming scenario, that is the stdevs and covs increase faster after 2050 than before (Figure 1).The next two sections explore the causes of the increased variability after global warming.

Strengthening Nutrient Control of Production Variability
In the Subarctic Atlantic, the summertime MLDs are shallow and the seasonal thermocline inhibits the supply of deep nutrients to the euphotic zone via vertical mixing.Yet, the supply of deep nutrients to the euphotic zone on an annual basis is large and dominated by wintertime entrainment, which occurs when light and hence productivity are low.Thus, on a climatological basis, the annual new production is approximately proportional to the physical supply of nitrate to the euphotic zone during the winter and the reduction of the nitrate in the euphotic zone between the late winter maximum and the late summer minimum (Figure 2a) (Whitt & Jansen, 2020;Williams et al., 2000).
Integrated over the Subarctic Atlantic, the March surface nitrate, annual new production, and growing-season chlorophyll anomalies are also approximately proportional at interannual and multidecadal timescales (Figure 2).The March surface nitrate variability can explain 40%-80% of the interannual new production and chlorophyll variance and 70%-100% of the multidecadal new production and chlorophyll variance (Figures 2b and 2e).In addition, there is approximately a mass balance between the March nitrate anomalies and the annual new production anomalies integrated over the top 50 m throughout the scenario, that is, the regressions on the interannual anomalies have slopes of about one (Figure 2c, blue line).The corresponding regression slope of the multidecadal anomalies is initially about 1.5 and declines to 1 with warming (Figure 2c, orange line).Thus, the March nitrate anomalies can only partially cause the correlated new production anomalies at multidecadal timescales early in the scenario.Yet, overall, the strong correlations and approximate mass balances support the nutrient-control hypothesis: that the physically mediated March nitrate anomalies are not merely correlated with the annual new production and chlorophyll anomalies but largely cause them.
As global warming progresses, nutrient control of integrated Subarctic Atlantic production and chlorophyll variability strengthens.The correlations of the regionally integrated March nitrate to production and chlorophyll  anomalies increase (Figures 2b and 2e).In addition, the local relationships between March nitrate and annual new production and growing season chlorophyll at scales down to 100 km strengthen and expand to include the regions of deepest March MLDs between 55°N and 65°N yielding more homogeneous dynamics across the Subarctic Atlantic (Figures S5-S7 in Supporting Information S1).As a result, the intensifying volatility of regionally integrated March nitrate (Figure 1h) ultimately causes elevated volatility in regionally integrated production (Figure 1e).However, because global warming also suppresses the sensitivity of chlorophyll to nitrate anomalies by a factor of 2 (Figure 2f), the stdevs of chlorophyll do not increase (Figure 1b) despite intensifying nitrate volatility and high correlations between nitrate and chlorophyll variability.The declining sensitivity of chlorophyll to nutrient anomalies is consistent with tightening top-down control of the phytoplankton biomass anomalies with global warming (Elsworth et al., 2022).The reduced variance of regionally integrated chlorophyll in the CESM (Figure 1b) thus reflects a balance between the increased variance of bottom-up drivers (Figures 1e  and 1h) and reduced sensitivity to those drivers (Figure 2f) in a warmer world.

Rising Sensitivity of Nutrients and Production to Physical Drivers
The changing physical climate is the ultimate cause of the changing biogeochemical variability.Thus, one might expect the increasing nutrient and production variability can be linked to increasing physical variability of some kind.Here, an alternative hypothesis is proposed: that the increasing volatility in the regionally-integrated nutrients and ultimately production in the Subarctic Atlantic is related to a rising sensitivity to variability of the AMOC and March MLD under global warming.Thus, the stdevs of nitrate and production can increase, while the stdevs of the AMOC and MLD decline as in Figure 1.

Interannual Timescales
Integrated over the Subarctic Atlantic, the March MLD anomalies are hypothesized to be the primary physical driver of the March nitrate anomalies in the euphotic zone on interannual timescales, because the March MLD is a measure of the annual entrainment of the higher-nitrate seasonal thermocline water into the surface mixed layer (Whitt, 2019;Whitt & Jansen, 2020;Williams et al., 2000).Complicating factors are that biogeochemical as well as advective processes also contribute to the temporal variability of the March nitrate concentration.In addition, there is spatial heterogeneity in the MLD and vertical nitrate gradient that rectifies on the sensitivity of the regionally averaged nitrate to the regionally averaged MLD in March.Nevertheless, the regionally averaged MLD anomalies can explain 60%-90% of the regionally averaged March nitrate variance in the euphotic zone on interannual timescales throughout the warming scenario (Figure 3a).And, interannual correlations between the March nitrate and MLD anomalies are found down to local scales of about 100 km throughout much of the Subarctic Atlantic (Figure S9 in Supporting Information S1).
In addition, global warming increases the sensitivity of the regionally averaged March nutrients to MLD anomalies.Specifically, the interannual correlations between the March nitrate and MLD anomalies increase with warming from about r 2 = 0.6 in 2020 to 0.8 in 2090 and the sensitivity doubles (Figures 3a and 3c).The increased sensitivity is qualitatively consistent with the forced reduction of the climatological March MLD (Figure 1j) and the strengthening of the nitrate deficit in the March mixed layer relative to the thermocline just below (not shown).
Consistent with strengthening nutrient control, the regionally averaged March MLD anomalies can explain a rising percentage of the regionally integrated production variance, from 10% to 20% in 2020 (notably less than 60% of the nitrate variance) to 70% in 2090 (Figure 3e; cf. Figure 3a).In addition, the sensitivity of interannual production anomalies to MLD anomalies more than doubles between 2020 and 2090 (Figure 3g), similar to the sensitivity of nitrate to March MLD variability (cf.Figures 3c and 3g).The increasing correlation and rising sensitivity are consistent with the expanding area of nutrient limitation and hence the area and magnitude of the positive correlation between the annual production and the March MLD and nitrate at local scales, including in the region of deepest March MLDs between 55°N and 65°N (Figures S5, S9, and S10 in Supporting Information S1).
Finally, there is little sensitivity of interannual production or nitrate anomalies to AMOC variability, consistent with the hypothesis of MLD control at interannual timescales (Figures 3a and 3e and Figures S11 and S12 in Supporting Information S1).

Multidecadal Timescales
In contrast to the higher-frequency interannual nitrate and production variability, the multidecadal variability is thought to be driven by both AMOC and MLD variability.Complicating factors include the changing biogeochemical rates and feedbacks in addition to rectifying effects of spatial heterogeneity on the regional integrals.Nevertheless, both the AMOC and regionally averaged March MLD variability are correlated with the regionally averaged March nitrate variability on multidecadal timescales (r 2 = 0.2 − 0.8; Figure 3b).And, multidecadal correlations of the March Global warming increases the sensitivity of regionally integrated nutrient anomalies to both AMOC and March MLD anomalies at multidecadal timescales.The sensitivity of the March nitrate in the top 50 m to AMOC anomalies is initially about 0.005-0.01mol N/m 2 per Sv (i.e., 0.1-0.2mmol/m 3 per Sv) from 2020 to 2050, after which it increases to about 0.02 mol N/m 2 per Sv in 2090 (i.e., 0.4 mmol/m 3 per Sv) (Figure 3d).The increased sensitivity of nitrate to MLD anomalies is qualitatively similar, from about 0.0005 to 0.0015 mol N/m 2 per meter between 2020 and 2090 (i.e., from 0.01 mmol/m 3 per meter in 2020 to 0.03 mmol/m 3 per meter in 2090).Interestingly, however, the correlations of both multidecadal AMOC and March MLD to nitrate increase non-monotonically with warming, from approximately r 2 = 0.5 at 2020 to r 2 = 0.2 at 2050 to r 2 = 0.8 at 2090 (Figure 3b).The increased correlations and sensitivities at 2090 (compared to 2020) occur in conjunction with stronger and more widespread correlations between the AMOC and the local March nitrate anomalies at scales down to 100 km (Figure S11 in Supporting Information S1).In addition, the regions of local correlation between the March MLD and nitrate anomalies shift spatially, from where the MLD is relatively modest in 2020 to where the MLD is deepest in 2090 (cf.Figures S9b and S9d in Supporting Information S1).The ultimately higher sensi tivities and correlations of nitrate to MLD or AMOC in 2090 are qualitatively consistent with a thinner and more nutrient-deplete seasonal thermocline relative to both the subtropical thermocline, where the upstream nutrients in the AMOC reside, and the thermocline just below the March MLD, the source of nutrients for entrainment anomalies (Whitt & Jansen, 2020).
Consistent with the nutrient control hypothesis, the changing correlations and sensitivities of the March nitrate to AMOC and the March MLD are approximately mirrored in the regionally averaged new production (cf.Figures 3b  and 3d to 3f and 3h).However, the local correlations to nitrate are systematically stronger and more widespread than the correlations to production (cf.Figures S11b and S11d to S12b and S12d in Supporting Information S1; and cf.Figures S9b and S9d to S10b and S10d in Supporting Information S1).
Future work is needed to fully understand and untangle the mechanisms of the changing multidecadal variability.The role of the MLD and AMOC are difficult to disentangle from a correlation analysis, because the March MLD variability is physically linked to AMOC variability on multidecadal timescales (Whitt & Jansen, 2020;Yeager et al., 2021).For example, Figure S13 in Supporting Information S1 shows that multidecadal variability of the regionally averaged March MLD is correlated with the AMOC variability and can explain the majority of the multidecadal AMOC variance after 2050.Moreover, the regions of strong local correlation between the March MLD and AMOC on multidecadal timescales are essentially the same regions where the March MLD is correlated with production and nutrient (cf.Figures S14b, S14d to S9b, S9d and S10b, S10d in Supporting Information S1).The more widespread local correlations of nutrient to AMOC than to the March MLD suggests AMOC may be the dominant driver of the multidecadal nutrient variability (c.f.Figures S11 and S12 to S9 and S10 in Supporting Information S1), but future work must separate the role of nutrient entrainment and meridional nutrient transport in multidecadal nutrient and production variability.Future work should also explain the non-monotonic responses to warming and low correlations around 2050.

Conclusions
This study adds to a growing body of research indicating that the Subarctic Atlantic Ocean biogeochemistry is especially and uniquely vulnerable to global warming.This study advances our understanding of this vulnerability by quantifying how regionally integrated phytoplankton production, biomass, and a controlling nutrient (nitrate) evolve in the high-emissions twenty-first-century RCP8.5 global warming scenario using a large ensemble of 32 simulations with the CESM.The ensemble enables relatively precise quantification of how the characteristics of internal interannual and multidecadal variability change under global warming.In fact, an ensemble of 32 runs is about the minimum size for detecting the statistically significant increases in the standard deviation of the annual new production between the 2020s and 2090s (p ≤ 0.05).However, a significant increase in the standard deviation of annual March nitrate anomalies (from the forced response/quadratic fit) in the 20-year windows centered at 2090 and 2020 is likely detectable from a single run (p ≈ 0.03).
The results confirm prior studies in that regionally integrated productivity, biomass, and nitrate decline significantly (by roughly half) in this scenario, in conjunction with similarly large fractional reductions in the March MLD and the strength of the AMOC, which together drive a substantial portion of the biogeochemical variability and change.
The new result is the finding that coefficients of variation for the regionally integrated production, biomass, and March nitrate in the Subarctic Atlantic all increase significantly (by roughly double between 2020 and 2090).Moreover, the standard deviations of annual new production and March nitrate in the euphotic zone both ultimately increase with warming, despite the major reductions in the corresponding means (Figure 1).An investigation of the causes of the interannual and multidecadal production variability supports a nutrient control hypothesis: that the regionally integrated production and biomass variance can largely be explained by physically-mediated variability in the regionally integrated March surface nutrients throughout the global warming scenario (Figure 2).An investigation of the physical drivers of the declining nutrient and hence production variability supports the rising-sensitivity hypothesis: that both the correlations and sensitivities of nutrients and production to physical forcing by MLD and AMOC variability increase with warming (Figure 3).Hence, surface nutrient and production variance ultimately increase, although the variance of both physical drivers declines in the warming scenario (Figure 1).
Since the response of Subarctic Atlantic biological production variability to global warming reflects an increasing sensitivity to declining physical variability, which is presumably reflected in its non-monotonic response to warming, future work should reassess the conclusion that global warming causes Subarctic Atlantic biological production variability to increase in other Earth system model large ensembles.In addition, future work should more fully investigate how ocean biogeochemical variability responds to global warming at the smaller spatial and temporal scales and in other regions than those considered here.

Data Availability Statement
The underlying CESM1-LE data sets that are the basis for this study are publicly available: https://www.cesm.ucar.edu/community-projects/lens/data-sets.The ensemble members used here are labeled 1-2, 9-24, 26-32, 34-35, and 101-105.For reference, the annual mean ecological nitrate tendency is labeled "J_NO3," the phytoplankton chlorophyll is a sum of contributions from all three phytoplankton groups (diatoms, small phytoplankton, and diazotrophs) "diatChl + spChl + diazChl," and the nitrate concentration is "NO3."The MLD in CESM is "HMXL."The AMOC is calculated from the monthly temperature ("TEMP") and salinity ("SALT") (converted to potential density using the model equation of state), the horizontal velocities ("UVEL" and "VVEL"), and nitrate concentrations.The horizontal advective nitrate fluxes are then projected onto density levels and then averaged in time and integrated along latitude lines.The visualization software and post-processed data are available at https://doi.org/10.5281/zenodo.8278032(Whitt, 2023).GlobColour data are available at http://www.globcolour.infoOceanColour-CCI data are available via http://www.oceancolour.org.

Figure 1 .
Figure 1.Measures of Subarctic Atlantic biogeochemistry (a-i) and physics (j-o) from the large ensemble of simulations with the Community Earth System Model version 1 in the RCP8.5 emissions scenario.The variables include (a-c) the May-August surface phytoplankton chlorophyll, (d-f) the annual new production integrated over the top 50 m, (g-i) the March nitrate averaged over the top 50 m, (j-l) the March mixed layer depth (MLD), and (m-o) the maximum of the Atlantic meridional overturning circulation (AMOC) streamfunction in density space at 45°N.The left columns include 20-year means centered every decade from 2020 to 2090 (•) and the range spanned by 32 ensemble members each year (shading).The middle and right columns show the standard deviations (stdevs, middle), the corresponding coefficients of variation (covs, right) normalized by the forced response, and their 95% confidence intervals (shading).The 1-year stdevs represent all years and ensemble members in eight overlapping 20-year windows from 2020 to 2090 (640 years per window), while the 20-year stdevs represent 32 20-year means from each ensemble member in each window (see Methods and Figures S1-S3 in Supporting Information S1).The α's in the middle column are confidence levels (in percentage terms: 100(1 − α)) that the stdevs are different in 2090 and 2020.

Figure 2 .
Figure 2. The anomalies in March nitrate integrated over the top 50 m are related to (a-c) the anomalies in the annual new production integrated over the top 50 m and (d-f) the anomalies in the May-August mean surface chlorophyll concentration (Chl).The scatter plots (a, c) show the relationships for all 95 years in all of the 32 ensemble members; the ensemble and 20-year means are overlaid (•).Then, the squared correlation coefficient (b, e) and the slope of the relationship (c, f) are derived from a linear regression in each 20 years window (centered every decade from 2020 to 2090).The shading in (c, f) represents the 95% confidence intervals.The regressions are calculated in two ways as described in the Methods: those labeled "1-year" represent interannual variability, whereas those labeled "20-year" represent multidecadal variability.The black lines in (c, f) indicate the slopes of the ensemble means through time (i.e., the slopes of the black circles in (a) and (d) respectively).

Figure 3 .
Figure 3. Relationships between the March nitrate (a-d) or annual new production (e-h) anomalies integrated over the top 50 m and the March mixed layer depth (MLD) anomalies (blue) or the Atlantic meridional overturning circulation (AMOC) anomalies at 45°N (orange).Squared correlation coefficients from linear regressions are shown in (a), (b) and (e), (f) for either the interannual "1-year" or the multidecadal "20-year" anomalies (see Methods).The corresponding slopes of the relationships are shown in (c), (d) and (g), (h) (the shading represents 95% confidence intervals).The thin red and blue horizontal lines in (c), (d) and (g), (h) indicate the slopes of the ensemble means through time, that is, the slopes that can be derived from the large color-filled black circles in Figures S8a and S8b in Supporting Information S1.