2.1 Time slice experiments

To develop a process-oriented means of examining the response of new production to idealized changes in climate, we adopt a time slice approach. Rather than running a full transient integration, we perturb the model's initial conditions and the surface forcing to simulate a period representative of the future climate state. We thus run separate integrations designed to be representative of early- and late-century climate conditions. For the former, we begin from a state initialized from observations and integrate the model with forcing representative of a statistically normal annual cycle, i.e., a normal year (Large and Yeager, 2004); we perform a 20-year physics-only spin-up, which is sufficient to reduce interannual drift in the physical state, and then use 10 further years, which include our biogeochemical model, as our early-century time slice. For our late-century time slice, we adjust the initial ocean state and atmospheric forcing variables using anomalies computed from the fully coupled CESM1 Large Ensemble (CESM-LE; Kay et al., 2015). The CESM-LE includes 40 members integrated from 1920–2100; we use anomalies computed from the ensemble mean, quantifying the difference in ocean state and atmospheric forcing variables between 2000 and 2100. We then integrate with the normal-year forcing plus monthly ensemble-mean atmospheric anomalies from the CESM-LE, again using 10 years as our time slice. The LE has been examined in comparison to both CMIP5 (Alexander et al., 2018) and observations (Deser et al., 2017), with results showing similarity to future sea surface temperature (SST) and observed climate variability, respectively. Using the century-scale mean changes from the CESM-LE allows us to represent the forced changes over the 21st century from the RCP8.5 scenario without the noise associated with natural interannual variability represented in any individual ensemble member.

Our resulting model runs have ocean temperature and salinity representative of early and late 21st century very similar to the LE. Drift in these values within the decadal runs is small compared to both the imposed climate perturbation changes between them and the typical interannual variability in a coupled model. The atmospheric surface state has the same sub-seasonal variability in each year and both epochs and no interannual variability. Physical fields of interest for the biological impacts of climate change include changes in available light, vertical stratification, and ocean currents. We discuss only these fields that can impact our idealized biogeochemical model, leaving out other fields that do not have direct impacts, such as temperature or pH. The changes in these fields are shown in Fig. 1, which may be summarized as follows: the Arctic receives more light (Fig. 1a), western boundary currents and the Antarctic Circumpolar Current speed up (Fig. 1b), equatorial upwelling spreads meridionally (Fig. 1c), winter mixing decreases in the Southern Ocean and shifts position in the North Atlantic (Fig. 1d), and near-surface stratification increases in most regions except the Arctic, with largest changes near the Equator (Fig. 1e). With this framework of the physical changes, we can consider their impacts on biological rates.

Figure 1Climate perturbation. (a) Change in maximum incoming short-wave radiation. (b) Change in annual-mean sea surface height (SSH). (c) Change in annual-mean vertical velocity at 100 m. (d) Change in annual maximum mixed-layer depth (MLD). (e) Change in annual-mean stratification at 100 m, $\mathrm{d}{\mathit{\sigma }}_{\mathit{\theta }}/\mathrm{d}z$, using 155 and 55 m potential density. All maps use a cylindrical equal-area projection.

2.2 Idealized tracers

2.2.2 2000s sensitivity

Our idealized tracers have five parameters: μ₀, k_N, α, σ, and w_s. In order to identify reasonable values, we perform a sensitivity analysis of the first three under the early-century climate. The decay rate and sinking rate of P are held fixed at w_s=0 m/d and $\mathit{\sigma }=\mathrm{1}/\mathrm{60}$ d in this analysis, as they do not affect N or the production rate and only modestly affect the horizontal spatial patterns of P: faster decay reduces global-mean P, while faster sinking moves P deeper in the water column, with neither changing the basin-scale or latitudinal structure of annual-mean surface P. The values of the parameters used are a factor of 4 increase and decrease from μ₀=0.5 mmol N/m³ d, k_N=1 mmol N/m³, and α=0.05 m²/W, giving three values for each. We ran 10-year simulations under 2000 conditions for each of the 27 cases created by combining the options and examined the zonal and annual-mean fields the final year; while this system is not in complete equilibrium after only 10 years, most drift occurs in years 1–2. We do not see substantial differences in our climate change results when using year 5 or year 10 of the time slices in our computations.

To examine sensitivity to parameters and choose reasonable cases for the climate change experiments, N is compared to near-surface World Ocean Atlas nitrate distributions (WOA NO₃, Garcia et al., 2013). P is compared to output from the CESM biogeochemistry model, BEC (Moore et al., 2013), which is run with the same model physics. Specifically, we compare near-surface P to total plankton, the sum of its three phytoplankton classes in carbon units, converted to nitrogen units using a stoichiometric ratio of 16 mmol N/117 mmol C. Although our particles represent both living matter and detritus, near the surface this P is most like newly produced plankton, which we expect to have the same spatial patterns as total phytoplankton. We choose this comparison to take advantage of the work done to fit BEC to observations for the same physics (see Fig. 4 of Moore et al., 2013) and due to the sparseness of phytoplankton concentration observations (e.g., Chl a data used in Boyce et al., 2010).

Overall, integrations with varied parameters for N mainly change its global-mean surface concentration and not the meridional structure (see Fig. 2a). This meridional structure matches that of surface global nitrate values qualitatively: high values in the subpolar regions, moderate near the Equator, and low in the subtropics. However, the Northern Hemisphere peak is shifted north, and the Southern Ocean values are low compared to observations, with the latter due mainly to our choice of a constant deep N concentration below that typical in the Southern Ocean. We will not be analyzing the Southern Ocean in detail, as we have not included iron limitation in our model, and this would be critical for the dynamics of production there. The annual and zonal mean nitrate from the BEC model are shown for reference (dashed line in Fig. 2a); this more detailed model also misses the location of the Northern Hemisphere peak, suggesting this may be due to the model circulation, but is much closer to observations in the Southern Ocean. Correlations between annual and zonal mean N and WOA NO₃ are between 0.81 and 0.88 for the full active depth of 1 km and between 0.65 and 0.93 for just the surface. Changes in parameters affect the mean surface concentration of N in a predictable way: increasing k_N increases the mean surface concentration, while increasing μ₀ and α decrease it. Changing μ₀ causes the largest change in magnitude and k_N the least.

This relationship between parameters and near-surface N is the result of a balance between physical supply of N from depth and its consumption by production at a rate set by the parameters. If production is slower than supply, N will accumulate, increasing the production rate and decreasing the vertical N gradient and therefore the supply rate until a near equilibrium is reached. Thus, we can predict this relationship between parameter choice and near-surface N using the production function, μ₀Q(N)L(I). The total production's dependence on N and I can be understood through the initial slopes of the production–nutrient and production–light curves, which are ${\mathit{\mu }}_{\mathrm{0}}/k_{\text{N}}$ and μ₀α: these describe how quickly production increases for an initial injection of each limitation; steeper initial slopes and faster initial production mean less accumulation of N before equilibration. Thus nutrient concentration will be higher for higher $k_{\text{N}}/{\mathit{\mu }}_{\mathrm{0}}$ and $\mathrm{1}/\left({\mathit{\mu }}_{\mathrm{0}}\mathit{\alpha }\right)$; these functions describe a two-dimensional space where increases along either axis decrease productivity and increase near-surface nutrient concentration.

We define a biological timescale as a single term to summarize these two axes. First, we recognize that $k_{\text{N}}/{\mathit{\mu }}_{\mathrm{0}}$ has units of days, so we normalize α, which has units W/m², by dividing by

$$\begin{array}{}\text{(5)}& {\mathit{\alpha }}_{\mathrm{0}}=\mathrm{1}({\text{m}}^{\mathrm{2}}/\text{W})\cdot (\mathrm{mmol}\mathrm{N}/{\mathrm{m}}^{\mathrm{3}})=\mathrm{1}\mathrm{mmol}\mathrm{N}/\mathrm{W}\mathrm{m}\end{array}$$

to reach the same units as $\mathrm{1}/k_{\text{N}}$, m³/mmol N, so that $\mathrm{1}/({\mathit{\mu }}_{\mathrm{0}}\mathit{\alpha }/{\mathit{\alpha }}_{\mathrm{0}})$ also has units of days. We then add the results. The choice of this normalization factor is specific to this model, where we are forming an equivalence between a nitrate-type nutrient and light. Our α₀ is an expression of how we stretch the light coordinate so that the initial slope of production with respect to I is in the same units as the initial slope of production with respect to N, suggesting α₀ as a ratio of $N/I$. Given the relative values of I and N, α₀≤1 mmol N/W m is likely to be the most fruitful.

Now we see that we have a biological timescale for new production based on the parameters chosen; we call this biological timescale τ_bio:

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{\text{bio}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\mu }}_{\mathrm{0}}}}\cdot (k_{\text{N}}+{\mathit{\alpha }}_{\mathrm{0}}/\mathit{\alpha }).\end{array}$$

The format of this timescale reinforces the qualitative relationship between surface N and changes in μ₀, α, and k_N. We can use τ_bio as a metric of how quickly production might consume a new supply of nutrients; if we compared it to a physical timescale for nutrient supply rate, τ_p, places where τ_bio is shorter are likely to be nutrient-depleted and places where τ_p is shorter are likely to be nutrient-replete. Across our 27 parameter cases, τ_bio ranges from about 2 d to 2 years. In Fig. 2b, τ_bio is compared to global-mean surface N for the 27 cases and correlates well, r=0.96. This correlation is best for this and similar values of α₀, e.g., 0.1 or 2 mmol N/W m, but is lower for, e.g., 0.01 or 100 mmol N/W m. This is consistent with our understanding of α₀ as a ratio of $N/I$.

Figure 2Annual and zonal mean of surface N (a) and P (c) for 27 parameter choice cases from the 10th year of simulation. Values for the two chosen cases are shown in blue and orange. Annual and zonal mean references are shown by the thick black lines: nitrate from the World Ocean Atlas observations and from the BEC model in the same physics (dashed); total plankton mass from the BEC model. (b) Case timescale values vs. global annual-mean N (black circles); the two chosen cases are marked with blue and orange stars.

Download

Meridional structures of P for these same integrations show more variation than N, which is clear in the more frequent intersections of the surface value curves and the variation in the latitudes of the maxima (Fig. 2c). The structure of the surface BEC phytoplankton has peaks near 60^∘ N, the Equator, 45^∘ S, and along the Antarctic coast. The three cases with $(\mathit{\alpha },{\mathit{\mu }}_{\mathrm{0}})=(\mathrm{0.0125},\mathrm{0.125})$ lack these peaks entirely, with a much smoother structure that we consider a poor match for observed behavior and will not be considered further. To determine the best of the remaining options, we measured the correlation coefficient between the annual meridional structure of P and BEC phytoplankton, which ranged between 0.18 and 0.68 at the surface. The top half of these cases, 12 of the remaining 24 (see Table 1), were examined by eye for a good qualitative match at the surface and will be used for further analysis. These cases have the correct general latitudinal structure of P and also have reasonable new (export) production rates between 5.5 and 7.6 PgC/yr (Fig. 3a), which is within the 5–11 PgC/yr rate in most literature (e.g., Cabré et al., 2015).

We will not be examining the explicit sinking of our particles. For those interested, consider sinking rates, w_s, consistent with estimates for detritus: single cells sink at about 0.1 m/d, and aggregates can reach over 100 m/d (Jackson and Burd, 1998). For decay rates, we suggest fixing these after choosing w_s and adjusting them to have most production sink below the depth threshold of choice. For instance, with w_s of 5 m/d, σ=1/yr has 95 % of the annual production sink below 100 m.

We choose two cases for detailed analysis which have P structures capturing different aspects of the BEC surface phytoplankton concentrations and quite different parameter and mean surface N values. The first case has surface P maxima near 45^∘ N and S and the Equator, nicely matching those at the Equator and 45^∘ S in BEC but missing the 60^∘ N and Antarctic-coast peaks. This case has a small maximum growth rate (μ=0.125 mmol N/m³ d), a moderate light sensitivity (α=0.05 m²/W), and a low nutrient threshold for growth (k_N=0.25 mmol N/m³). We consider this analogous to a small phytoplankton type or a phytoplankton community that has adapted to oligotrophic conditions, with plenty of light but low available nutrients. This case has a τ_bio of 162 d, corresponding to moderate mean surface N, slightly lower than WOA nitrate at the surface (see Appendix A for more details). We call this the “slow” case.

The second case has surface P maxima near 60^∘ N, the Equator, and the Antarctic coast, similar to those at 60^∘ N and the Equator in BEC and capturing the increase toward the Antarctic coast but missing the 45^∘ S maximum. In contrast to the first, this case has a fast maximum growth rate (μ₀=2 mmol N/m³ d), a low light sensitivity (α=0.2 m²/W), and a moderate nutrient threshold for growth (k_N=1 mmol N/m³). We consider this analogous to a large phytoplankton type or a phytoplankton community that has adapted to higher-latitude conditions, with large seasonal cycles of light and nutrient availability. This case has a τ_bio of 3 d, corresponding to a very low mean surface N, which is noticeably below observed nitrate concentrations. We choose this case for large contrast with the first and call it the “fast” case. In the results section, we will focus on the changes in production under global warming for these two parameter cases of our idealized biogeochemical model and use the larger set of 12 cases for context.

Table 1Set of parameter values for 12 cases that compare best to WOA nitrate and CESM total phytoplankton. Each column represents three cases with $k_{\text{N}}\in (\mathrm{0.25},\mathrm{1},\mathrm{4})$. The range of τ_bio covers those for all three k_N.

Download Print Version | Download XLSX

2.3 Light and nutrient controls on production

Under global warming, changes in light and nutrient availability will vary both spatially and temporally. This section describes our method to untangle the effects of these two components on changes in new production. The form of production, R=μ₀QL, allows for a decomposition and attribution of changes in productivity to changes in nutrients and/or light. As μ₀, the maximum growth rate, does not vary in space or time, we can examine simply the nutrient availability, Q, and the light availability, L, both of which are nondimensional and have values between 0 and 1, as does their product. Using model output of the monthly-mean N and R, we compute the monthly-mean $Q=N/(k_{\text{N}}+N)$ for each grid cell and the monthly-mean $L=R/{\mathit{\mu }}_{\mathrm{0}}Q$. The reason for computing L in this way is that the mixed-layer depth, and thus L, can change rapidly. The production rate, R, and nutrient, N, are averaged online in our model, and thus computing L from them allows us to have a time-averaged L that would not be possible to compute from the time-averaged mixed-layer depth and incoming radiation. The change in QL can be decomposed as

where Δ indicates the change between the warmer, perturbed climate (2100) monthly values and those in the 2000s climate, as written out using subscripts of 2100 and 2000, and all terms can be averaged in space or time as needed. These difference terms are also nondimensional and have values between −1 and 1. Due to the way we compute L from R and Q, there is no residual term – the above equation holds to numerical precision. We will use this decomposition to analyze both the spatial patterns of annual-mean production and seasonal cycles of production. When one of the three right-hand terms is most similar in size to ΔQL and highly correlated in space or time, we will consider that term the main driver of changes in production.

For spatial correlations, or pattern correlations, we assume a decorrelation length scale of 10 grid cells, which is 6–9^∘ latitude or longitude, depending on location. From examining the autocorrelation of several terms in Eq. (7), we find that this is a slight overestimate of the decorrelation length scale in some cases. The degrees of freedom in a pattern correlation over the global ocean are then 825, and pattern correlation coefficients larger than r=0.09 are significant.

Figure 3(a) Global annual new production integrated over the top 100 m for 12 cases for 2000 (black circles) and 2100 (red diamonds) vs. τ_bio (nonlinear axis for readability). (b) Percent change from panel (a), black stars; fast and slow cases circled in orange and blue, respectively. (c) Global mean N in the top 100 m vs. τ_bio. (d) Absolute change from panel (c).

Download

3.1 Global statistics and dominant spatial patterns

For our 12 cases, global annual new production in the top 100 m is 5.3–7.5 PgC in the early-century scenario and decreases to 4.8–6.2 PgC in the late-century scenario; production is higher for shorter τ_bio in both epochs (Fig. 3a, b). These patterns hold for our two exemplar cases, with the slow case having lower new production and a smaller decrease than the fast case. Our 12-case range of decreases in production, 9.5 %–19.5 %, is similar to the 7 %–18 % range of export decreases in CMIP5 (Fu et al., 2016) but larger decreases than most net primary production changes in CMIP5, −2 % to −16 % (Fu et al., 2016), or CMIP6, +6 % to −12 % (Kwiatkowski et al., 2020). A large portion of this variability in global new production is related to the Southern Ocean, which is the basin with largest production and which our model does not represent well. Without the Southern Ocean the reductions in global annual production are 8.5 %–11 % – a range which is smaller than the range of export decreases in CMIP5 and within the range of net primary production changes in both CMIP5 and CMIP6.

Decreases in new production are partially determined by decreases in near-surface nutrient concentration. In all 12 cases, global-mean N concentration in the top 100 m decreases by 15 %–22 % in the late-century conditions (Fig. 3c, d). Both initial concentrations and absolute reductions are smaller for shorter τ_bio, but the small range of reduction percentages (15 %–22 %) highlight that the changes are somewhat insensitive to the varied light- and nutrient-limitation choices. These absolute reductions are 0.04–0.66 – a range which is slightly smaller than the reductions of 0.66±0.49 for CMIP5 RCP8.5 and 1.06±0.45 for CMIP6 SSP5-8.5 (Kwiatkowski et al., 2020). Despite the small absolute reductions in near-surface N concentration, the decrease in global new production is larger for fast cases with the shorter case timescale (Fig. 3b); these shorter τ_bio cases have lower near-surface nutrient concentrations in the early-century conditions and higher slopes in their production–nutrient curve, making them more sensitive to changes.

The spatial patterns of annual production under early-21st-century conditions demonstrate one effect of τ_bio (Fig. 4a–b). The slow case has lower maximal values, as expected given its lower maximal growth rate, and smoother variations in space. The fast case has higher maximal values of new production, again set by its maximal growth rate, and sharper variations in space, due to its greater need for nutrients in order to grow: production is more limited to locations with a fast physical nutrient supply. Despite these differences, both cases match the general meridional patterns of observed phytoplankton (see, e.g., Dasgupta et al., 2009 for observations), with higher values in the subpolar and equatorial regions and lower values in tropical–subtropical regions.

Figure 4(a, b) 2000s annual production integrated over the top 100 m. (c, d) Percent change with climate perturbation. (a, c) Slow case. (b, d) Fast case. The magenta contour indicates the subtropical South Pacific, cyan the Porcupine Abyssal Plain, and the black line the southern edge of the Arctic.

The spatial patterns of change in production under the warmer climate are qualitatively similar for both cases, with broad moderate reductions over most of the oceans, including the Indian Ocean, the South Pacific, and the North Atlantic (Fig. 4c–d; the pattern correlation is r=0.26; see Table 2 for all pattern correlations). Both cases also contain common areas with increased production in the 2100s, including most notably the edges of the equatorial Pacific, where upwelling has spread latitudinally, and the southern edge of the subtropical North Pacific Gyre. Despite broad similarity in the productivity response between the fast and slow cases, there are many small regions where the sign of the climate change response differs, perhaps most notably the Arctic, the central equatorial Pacific, and the Bay of Bengal. These are regions with qualitative as well as quantitative sensitivity of the climate perturbation reaction to the light and nutrient functions of uptake.

While we do not discuss spatial patterns in detail for the larger set of parameter cases, we do compute production by basin (Appendix, Fig. B1). For every basin but the Arctic, production decreases for a warmer climate in all 12 cases, indicating that the qualitative results in most basins are not sensitive to the specifics of light and nutrient limitations. In the Arctic, which has the least total production, the four cases with lowest τ_bio have decreases in new production, while the remaining eight have increases. We will discuss this sensitive region further in Sect. 3.5.

As global reductions in production were matched by reduced nutrient concentrations, so too can spatial patterns in production be linked to nutrient concentrations and the vertical nutrient flux. Although the nutrient flux and new production do not necessarily align in space and time, we empirically find pattern correlations between the annual-mean production and the nutrient flux at 100 m as well as the mean nutrient in the top 100 m (Fig. 5). The vertical flux is somewhat more noisy; hence the correlations of the absolute changes between present and future are only r=0.51 and r=0.14 for the fast and slow cases, respectively (see Table 2 for all pattern correlations). The mean nutrient concentration in the top 100 m reflects the integrated effects of the flux, the stronger vertical mixing in the mixed layer, and productivity; hence the correlations are somewhat stronger. We note that the variations in nutrient concentrations within the mixed layer are typically small – less than half the mean concentration and much less than k_N; these variations lead to slightly lower production rates than those that would occur if the nutrient concentrations were constant from the surface to the mixed-layer depth. In the present day, the pattern correlations between annual nutrient concentration in the top 100 m and production integrated to 100 m are r=0.56 and r=0.26 for the fast and slow biogeochemistry; for the future they are r=0.58 and r=0.33, and for the percent changes they are r=0.88 and r=0.48, respectively. Across our 12 cases, the lower τ_bio cases, which are more nutrient limited, have stronger correlations between nutrient concentration in the top 100 m and productivity. These correlations are always stronger in the warmer climate, when all cases are more nutrient limited.

Figure 5(a, b) Nutrient concentration, mmol N/m³, averaged over a year and the top 100 m. (c, d) Percent change of this field from 2000s to 2100s climate. (a, c) Slow case. (b, d) Fast case.

While we have connected the production and its changes to nutrient concentrations, light effects are also important. Productivity is a product of functions representing the sensitivity of productivity to light L(I) and to nutrient Q(N) availability. As described in the methods, this allows a decomposition of the changes in production into those caused by changes in nutrient availability, LΔQ; those caused by changes in light availability, QΔL; and the covariance of the two, ΔQΔL. Examining the spatial patterns of QL, its change (Δ(QL)), and the components of that change (Fig. 6) provides more details on the controls of reduced production beyond the vertical nutrient supply and surface concentration. From the spatial fields of these terms, all annual means averaged over the top 100 m, QL, and Δ(QL) show the same spatial patterns as production and its changes, respectively, as expected by definition. However, Δ(QL) looks somewhat different from the percent changes in production shown before (Fig. 4). From the components of the change, it is clear that LΔQ is the main control on production: r=0.74 and r=0.64 for the fast and slow cases, which is consistent with our discussion of nutrient supply and concentration (Fig. 6b, c, g, h). The change in light availability is a smaller contributor: r=0.19 and r=0.43 for the fast and slow cases (Fig. 6d, i). Finally, the covariance is anticorrelated to the total change: $r=-\mathrm{0.054}$ and $r=-\mathrm{0.025}$ for the fast and slow cases (Fig. 6e, j); these correlation coefficients are insignificant. This component is largest near the Equator, where it offsets increases in nutrients due to broader upwelling in the warmer climate. These pattern correlations are qualitatively similar without the Southern Ocean, which we do not represent well with this model.

Figure 6QL for 2000s climate (a, f), its change (b, g), and its components, all averaged over 1 year and the top 100 m and then normalized by the maximum of QL (nondimensional) in the 2000s. (c, h) LΔQ, (d, i) QΔL, and (e, j) ΔQΔL. (a–e) Slow case, normalized by 0.055. (f–j) Fast case, normalized by 0.14. Minima and maxima for each plot are given following the panel letter.

From these analyses, we see that qualitatively the spatial patterns of the response to the climate perturbation are similar across the biogeochemical cases under consideration, likely pointing to strong controls by the underlying physics via the nutrient supply. Quantitative differences exist; the spatial differences point to locations where results are sensitive to model parameters, while the global differences are consistent with the chosen parameters in that the faster cases, which have higher production and lower near-surface nutrients, have new production more correlated with the reductions of near-surface nutrients and their vertical supply.

Table 2Pattern correlations for annual-mean fields. Bold values indicate significance. All production, production component, and nutrient concentration (N) fields are the mean over the top 100 m, and all fluxes are at 100 m depth.

Download Print Version | Download XLSX

3.2 Global seasonal cycle

Over much of the ocean, the seasonal cycle is the dominant mode of productivity variability due to seasonal variations in light and nutrients. We are interested in how the seasonal cycle will respond to our climate perturbation and whether this is sensitive to the case parameters.

We define the global seasonal cycle as monthly global spatial averages with a 6-month offset between the Northern Hemisphere and Southern Hemisphere so that the boreal and austral seasonal cycles are in phase. Differences in the seasonal cycle of new production between the fast and slow cases are large and are a reflection of the different τ_bio of the cases (Fig. 7b). Seasonal changes in nutrient and light availability have physical timescales of weeks to months, which are between the τ_bio of 3 d for the fast case and 162 d for the slow case. Both cases have high upward nutrient flux in the winter and early spring (Fig. 7a) in conjunction with deep and dense winter mixed layers. However, the strong nutrient flux occurs over a shorter period of time in the slow case as the physical supply overwhelms the slow ecosystem's ability to consume and thereby reduces the flux. Productivity in the fast case quickly consumes the nutrients as they are supplied in the late winter and spring, whereas productivity in the slow case spreads out the nutrient consumption, achieving its maximum in midsummer when the light is most plentiful (Fig. 7b).

Despite the qualitative differences between the present-day seasonal cycles, climate change has a qualitatively similar impact on both the fast and slow cases (Fig. 7a–c). The peak wintertime nutrient flux at 100 m (Fig. 7a) is reduced by 22 % and 21 % in the fast and slow cases, respectively. This reduction is consistent with, but may not be entirely caused by, reduced winter mixed-layer depths restricting the entrainment of nutrient-rich water from the thermocline (Fig. 7f). Production is reduced during all months but particularly in the latter half of the respective growing seasons. Thus, the growing seasons during the 2100s are modestly shorter with earlier peaks in both cases. It is true across all 12 cases that the largest reductions in new production are after the peak new production. Thus, we can conclude that while the seasonal cycle can be very different across cases, for this model the global reduced production in a warmer climate follows a consistent pattern of a shortened growing season.

The causes of reduced production can be identified through the seasonal cycles of the components of ΔQL, all of which are small fractions of the maximum early-century QL values. In both slow and fast cases, winter and spring increases in light availability due to shallower MLD are offset by reductions in the nutrient availability and light–nutrient covariance, leading to reduced production (Fig. 7d, e). The minima of ΔQL are in the second half of the growing season in both cases, driving the shortened season. The slow case has negative values of all components at this time, with nutrient availability being the most negative (Fig. 7d). By contrast, the fast case's minimum of ΔQL occurs while QΔL is positive and reduced production is driven by both reduced nutrient availability and the light–nutrient covariance (Fig. 7e). This large ΔQΔL in the fast case during the growing season is obscured in the annual mean (Fig. 6), where this factor is small over most of the ocean and insignificantly correlated to the total change. In both our cases, reduced nutrient availability is a major contributor to the shortened growing season, indicating a consistent mechanism for reduced total production in response to the climate perturbation.

Figure 7Globally integrated supply of nutrients (a) and production of particulates (b) for the top 100 m. Both cases are shown for the 2000s and 2100s; the seasonal cycle is the mean of simulation years 8–10. Production controls (nondimensional) for the slow case (c) normalized by a global maxQL(2000) value of 0.055 and the fast case (d) normalized by a global maxQL(2000) value of 0.14. (e) Mean mixed-layer depths.

Download

3.3 Regional case studies

We next examine three regions in detail: the downwelling South Pacific, the Arctic, and the Porcupine Abyssal Plain in the North Atlantic. These three regions exemplify the spatially varied processes driving the climate response and the sensitivity of that climate response across our model cases (Figs. 4, 6). We emphasize that none of the model cases closely match observations in any particular region (see Appendix A for a comparison of N and observed nitrate) and that these examples are chosen for their qualitative differences. To contextualize these regions, we note their physical biomes; biome analysis has been used previously to identify ocean regions with similar biogeochemically relevant physical characteristics across basins (Sarmiento et al., 2004; Cabré et al., 2015). A further discussion of biome analysis across the broader 12 parameter cases is available in Appendix B.

The downwelling South Pacific (red outline in Fig. 4), part of the permanently stratified subtropical biome, is a region with small seasonal cycles, small climate-induced changes in physics, and broad losses in production across model cases. These losses of production are driven by the same processes in the fast and slow cases, suggesting an insensitivity to parameter choices. By contrast, the Arctic (north of 66.5^∘ N, the black line in Fig. 4), in the ice and marginal ice biomes, has large climate-induced physical changes, mainly in ice coverage. As noted earlier, the climate response in the Arctic is sensitive to parameter choice; this region is the most sensitive of those we discuss, with changes in annual new production having opposite signs in our two cases. The Porcupine Abyssal Plain (blue outline in Fig. 4), in the seasonally stratified subtropical biome, has deep winter mixing in the present climate and large decreases in mixed-layer depths with the climate perturbation. This region is an example where the main driver (light or nutrients) of the climate-change-induced reductions to productivity changes across parameter space. These three regions provide illustrative examples demonstrating how physical climate changes create a range of possible responses in new production, as mediated by the parameters that modify the sensitivity of productivity to light and nutrient availability, and thus how the sensitivity of projected changes in new production varies across the world oceans.

3.4 South Pacific

The downwelling South Pacific is defined by annual-mean downward vertical velocity at 100 m depth in the South Pacific between 10 and 35^∘ S. Maximum winter mixed-layer depths reach only 120 m, putting this region mainly in the permanently stratified subtropical biome. Here, all biogeochemical parameter cases are qualitatively similar in both their current-climate seasonal behavior and their response to the climate perturbation, which we show is driven by reduced vertical nutrient supply.

The seasonal cycles in this region in both early- and late-century climates have high values of both upward nutrient flux and new production in the austral winter and spring (Fig. 8a, b). Peak upward nutrient flux is in August in both cases in the early- and late-century conditions. In the fast case, peak new production is in the same month, indicating a fast response to the supply of nutrients. In the slow case, there is a 1-month delay, indicating a slower response. These peak times do not change between early- and late-century conditions; only the overall magnitude of rates decreases: in the slow case, annual production decreases by 25 % from 0.11 to 0.086 PgC/yr, while in the fast case, annual production decreases by 39 % from 0.045 to 0.027 PgC/yr. This is consistent with the global production being more reduced in the fast case than in the slow case. This pattern holds across the 12 cases, with all having peak production in August through October, no change in the peak timing with climate, and larger decreases in production for shorter τ_bio.

The decrease in new production is due to the reduced availability of nutrients. The total change in production, Δ(QL), is nearly identical to the changes due to nutrient availability, LΔQ, with rms differences of $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ for the fast and slow cases (Fig. 8c); this near match holds across all 12 parameter cases. This is consistent with Cabré et al. (2015), who saw decreases in production in low-latitude nutrient-limited biomes across models. Our addition to that analysis is the direct quantitative connection to reduced nutrient availability that is not straightforward to compute for more complex biogeochemical models.

The reduced nutrient availability is due to reduced upward fluxes of nutrients, which we can exactly identify (Fig. 9a, b). These fluxes are comprised of downward advective fluxes, upward fluxes from parameterized vertical mixing (Large et al., 1997), and upward fluxes from vertical effects of parameterized along-isopycnal mixing (Redi, 1982; Gent et al., 1995). Both components of parameterized mixing show reduced annual fluxes in the warmer climate at 100 m. Neither mixed-layer depths nor effective vertical diffusivities change substantially – less than 10 % for all MLD and less than 10 % across $\mathrm{2}/\mathrm{3}$ of the grid points in this region for diffusivity. Thus, reduced $\mathrm{d}N/\mathrm{d}z$ is the driver of the decreased upward nutrient flux.

The changes in spatial-mean nutrient profiles (Fig. 9c) show a consistent pattern across the fast and slow cases: increased nutrient at 350–650 m depth but decreases from 350 m to the surface, including decreases near 100 m. These mid-depth nutrient changes are due to physical changes in circulation and/or mixing in the main thermocline. For diffusive transport, the decreases in N concentration above 350 m mean that the vertical gradient of N near 100 m, where we measure the fluxes, decreases. Thus, the decreased production in the downwelling South Pacific is being driven by decreased nutrients below the deepest winter mixed layers (about 120 m), through decreased $\mathrm{d}N/\mathrm{d}z$ causing decreased diffusive fluxes, decreased near-surface nutrients, and thereby decreased nutrient availability, Q. From these two examples it seems likely that the causes of reduced nutrient availability in the warmer climate, which drives reduced production, are consistent across parameter cases.

Figure 8Subtropical South Pacific nutrient supply (a) and production of particles (b). Production controls (c) normalized by maxQL(2000) values in the region: 0.032 for the slow case and 0.013 for the fast case.

Download

Figure 9(a–b) Annual-mean upward N flux components, gC/m² yr, at 100 m depth, averaged over the subtropical South Pacific, for the 2000s and 2100s climate and their difference. Parameterized vertical mixing (“vert mix”), vertical component of parameterized isopycnal mixing (“iso mix”), and advection (“adv”) from model diagnostics. (a) Slow case. (b) Fast case. (c) Profiles of change (2100–2000) in annual-mean subtropical South Pacific nutrients.

Download

3.5 Arctic

The Arctic region is defined as being within the Arctic circle (north of 66.5 ^∘N) or equivalently having at least 1 d per year with no incoming solar radiation. This region has the largest sensitivity of changes in new production to model light and nutrient limitation, with the range across 12 parameter cases being a decrease of 17 % up to an increase of 52 %. Here, our two example biogeochemistry cases have notably different responses to climate change, with the slow case having a 46 % increase in annual production and the fast case having a 16 % decrease.

Seasonal cycles in the 2000s are similar for the fast and slow cases. In this region, that is an upward flux of nutrients year-round and high new production rates in the summer, with a peak in May for the fast case and September for the slow case (Fig. 10a, b). Under late-century conditions, the peak production of the fast case is in the same month, with a slight increase in production earlier in the year and reductions later in the year which lead to a total decrease; in contrast, the peak production of the slow case is much higher and is earlier – now in July. These regional-mean production shifts are consistent with the different responses seen in the annual production maps (Fig. 4), where the slow case has large increases while the fast case has moderate changes of both signs.

In the Arctic, increases in incoming light, with the seasonal maximum of the regional mean more than doubling (91 to 202 W/m²), are due to reduced sea ice (Fig. 1a) and have a large impact on changes in production. The mean increase in light availability in summer, quantified by QΔL, its effect on changes in production, is nearly as large as the maximum QL in the 2000s for both cases (Fig. 10c, d). In the slow case, this increase in light allows for a large increase in production in the first half of the growing season, until nutrients are slightly depleted by that production, reducing production through lower nutrient availability in the second half (negative LΔQ and ΔQΔL). This pattern is qualitatively consistent across slow cases with long τ_bio. In contrast, the increase in light availability in the fast case is offset by an associated reduction in nutrient availability, such that $\mathrm{\Delta }Q\mathrm{\Delta }L\approx -Q\mathrm{\Delta }L$. Increased light in the spring leads to immediate increases in production (April), which uses enough nutrients to cause a dip in nutrient availability (May) before peak light availability (July), leading to decreased production in the summer and fall. This pattern is consistent across the four cases with shortest τ_bio, which have decreased annual production in the warmer climate (Fig. B1). Thus, in the Arctic, the increases in light availability consistently increase production in spring, leading at some point to a biologically driven decrease in near-surface nutrients. The resultant changes in production may take either sign and depend on the speed at which nutrients are removed, consistent with our τ_bio.

The Arctic is a region where changes in new production under global warming is highly sensitive to the formulation of model production, both for our model and more complex ones. Vancoppenolle et al. (2013) found that CMIP5 projections of Arctic NPP were dependent on whether the Arctic reached a nutrient-limited state post-sea-ice loss, with those models that reach a nutrient-limited state having reduced production and others having increased production in the warmer climate. Our analysis provides a mechanistic hypothesis for these differences: shorter τ_bio models have lower near-surface nutrients in early-century conditions, and therefore increases in production with increased light will more quickly lead to a nutrient-limited state.

Figure 10Arctic nutrient supply (a) and production (b). Production controls on the slow (c) and fast (d) cases, normalized by the maximum values of arctic maxQL(2000): 0.021 for the slow case and 0.028 for the fast case.

Download

3.6 North Atlantic

The Porcupine Abyssal Plain in the northeast North Atlantic is defined here to be 40–52^∘ N and 27–11^∘ W (cyan box in Fig. 4); this includes the location of the Porcupine Abyssal Plain Sustained Observatory. This region is characterized by deep winter mixed layers and net downwelling. Under a warmer climate, the winter mixed-layer depths are reduced; the spatial mean of the maximum decreases from 280 to 229 m, and the absolute maximum MLD is reduced from 669 to 590 m (Fig. 11d). All 12 cases have decreased new production in the future, 9.55 % to 26.75 %, which is a smaller range than our other two regions but still larger than the global ocean; our two example cases have a decrease of 12 % for the slow case and 21 % for the fast case. In this region, as for the Arctic, the mechanisms driving the changes vary.

In this region, our two model cases show qualitatively different seasonal cycles with production mainly in the early spring for the fast case and spread over the summer for the slow case (Fig. 11b) similar to their respective global cycles. Both cases have little qualitative change between the early and late 21st century. Most of the nutrient supply is in the winter for both cases (Fig. 11a), mainly due to wintertime entrainment associated with the mixed layer.

In the warmer climate, changes in production, ΔQL, depend on the relative impacts of increased light availability, QΔL, in March and April and decreased nutrient availability, LΔQ, in all months; the signs of these drivers are consistent across cases, but the signs of the change in production are not. In the slow case, the changes in production with climate are an increase in spring and a decrease in summer and fall for a total decrease, closely following changes in light availability, Δ(QL)≈QΔL (Fig. 11c). The increased light in March and April corresponds to shallower mean and maximum mixed-layer depths, respectively. Light availability decreases in the summer, due to increased mixed-layer depths in May and June. In the summer and fall, a higher portion of the reduction in production is due to reduced nutrient availability from both the lower winter peak in nutrient supply and the increased use during the spring production.

In the fast case, the same physical changes of the climate perturbation result in production decreases due to reduced nutrient availability in all months, Δ(QL)≈LΔQ (Fig. 11c). The winter increase in light availability has no noticeable impact, with substantially lower winter–spring production co-occurring with the lower nutrient flux (Fig. 11a) and shallower monthly-mean mixed-layer depth (Fig. 11d). This region's changes in new production are sensitive to parameter choices. Although total new production is reduced in both cases, reduced winter mixed-layer depths can act to either increase or decrease spring production, depending on whether that production is more sensitive to light, as in the slower cases with longer τ_bio, or nutrient availability, as in the faster cases with shorter τ_bio. These differences highlight the usefulness of this idealized model, which allows us to diagnose these drivers.

Figure 11Porcupine Abyssal Plain nutrient supply (a) and production of particles (b). Production controls (c), normalized by maxQL(2000) values for the 2000s: 0.022 for the slow case and 0.015 for the fast case. The discrepancy for the slow case is mainly due to LΔQ. Range and mean of monthly-maximum MLD (d).

Download