A latitudinally banded phytoplankton response to 21 st century climate change in the Southern Ocean across the CMIP 5 model suite

Changes in Southern Ocean (SO) phytoplankton distributions with future warming have the potential to significantly alter nutrient and carbon cycles as well as higher trophic level productivity both locally and throughout the global ocean. Here we investigate the response of SO phytoplankton productivity and biomass to 21st century climate change across the CMIP5 Earth System Model suite. The models predict a zonally banded pattern of phytoplankton abundance and production changes within four regions: the subtropical (∼ 30 to 40 S), transitional (∼ 40 to 50 S), subpolar (∼ 50 to 65 S) and Antarctic (south of∼ 65 S) bands. We find that shifts in bottom-up variables (nitrate, iron and light availability) drive changes in phytoplankton abundance and production on not only interannual, but also decadal and 100-year timescales – the timescales most relevant to climate change. Spatial patterns in the modelled mechanisms driving these biomass trends qualitatively agree with recent observations, though longer-term records are needed to separate the effects of climate change from those of interannual variability. Because much past observational work has focused on understanding the effects of the Southern Annular Mode (SAM) on biology, future work should attempt to quantify the precise influence of an increasingly positive SAM on SO biology within the CMIP5 models. Continued long-term in situ and satellite measurements of SO biology are clearly needed to confirm model findings.


Introduction
The photosynthetic activity of marine phytoplankton provides the ultimate source of food for virtually all marine biota, including organisms of vast commercial value.This phytoplanktonic activity also drives the biological pump, the process by which surface carbon dioxide and nutrients are drawn down via photosynthesis with subsequent sinking of organic matter to the deep ocean that effectively removes carbon from the atmosphere for centuries to millennia (Eppley and Peterson, 1979;Heinze et al., 1991).The warming trend recorded in the global surface ocean since the mid-20th century is projected to continue in the 21st century (Stocker et al., 2013) and can impact phytoplankton activity both directly via the physiological effect of temperature on growth rate and/or indirectly by altering key environmental factors such as nutrient and light availability (e.g.Marinov et al., 2010).The responses of phytoplankton communities to climate change may have profound ecological and biogeochemical repercussions with potential feedbacks on climate, the net sign and magnitude of which are still largely uncertain.Documenting and understanding these responses is one of the main goals of global change science today (Falkowski et al., 2000;Geider et al., 2001).
As a major region of deep, intermediate and mode water formation, the Southern Ocean (SO) is one of the few places on Earth where there is direct communication between the atmosphere and the deep ocean.Because of this, the SO plays a critical role in the global climate system via its significant impacts on the global heat and carbon budgets.Additionally, intermediate and mode waters formed here allow for large advective transfers of macronutrients such as nitrate, phos-phate and silicate from the SO to the low-latitude oceans, indirectly accounting for up to 75 % of phytoplankton production north of 30 • S (Sarmiento et al., 2004a;Marinov et al., 2006).Thus, potential changes in SO productivity can affect not only local nutrient and carbon cycles, but may also drastically alter nutrient and carbon cycles as well as phytoplankton distributions throughout the global ocean.
Much of the SO is a so-called HNLC (high-nutrient, lowchlorophyll) region, where chlorophyll concentrations (and implicitly phytoplankton biomass and production) are relatively low, in spite of a large upwelled supply of macronutrients (e.g.Martin et al., 1990;Cullen, 1991;Pitchford and Brindley, 1999).Here insufficient light availability may help explain why biological productivity is not as high as it could be.Because light is a potentially stronger limiting factor than macronutrient supply for photosynthesis here, warming is generally postulated to be advantageous for algal communities within these regions because shallower mixed layer depths (MLDs) (due to enhanced stratification and increased freshwater influx with future warming) are expected to increase light availability to phytoplankton and prolong the growing season (Bopp et al., 2001;Le Quéré et al., 2005;Doney, 2006).Warming may also directly enhance productivity by alleviating growth rate limitations due to low temperatures (Steinacher et al., 2010).If this line of reasoning holds, we should observe an increase in phytoplankton biomass and chlorophyll concentrations in the high-latitude SO with future warming.A further complicating factor, however, is that SO phytoplankton are also limited by iron and silicate, such that they can be light-iron-silicate (or any combination of the three) co-limited (C.M. Moore et al., 2013).Thus, changes in any of these factors will affect phytoplankton productivity and biomass within the SO.Because of the complicated multifactorial nature of the problem, a synergy of observations and models is needed to understand the driving mechanisms of projected changes in SO phytoplankton distributions.
Recent studies have suggested that SO phytoplankton biomass and productivity will change in response to rising atmospheric CO 2 concentrations, but the direction, significance, and causes of these changes are still under debate (Bopp et al., 2001(Bopp et al., , 2005(Bopp et al., , 2013;;Schmittner et al., 2008;Steinacher et al., 2010;Wang and Moore, 2012;Marinov et al., 2013;Cabré et al., 2014;Laufkötter et al., 2015).Here we use the newest generation of fully coupled CMIP5 (Coupled Model Intercomparison Project 5) Earth System Models to systematically study the response of SO phytoplankton to 21st century climate change, assuming the rcp8.5 emissions scenario.To this end, we borrow some statistical methods developed in Cabré et al. (2014) (namely, the model weighting scheme and the bootstrap technique, both described in Section 2 below) to conduct our work.All 16 of the CMIP5 models that incorporate ecological subroutines and provide their output on the CMIP5 portal are included in our study.We also summarize and review past field studies of SO phy-toplankton to see what has already been observed and to understand where there may be disagreement over mechanisms and/or recent directions of changes between the models and field data.We find that over the next 100 years, the CMIP5 models predict a zonally banded pattern of SO phytoplankton abundance and productivity changes driven by shifts in light, nitrate and iron availability with future warming.We show that the SO south of ∼ 30 • S can be separated into four zonally defined biomes: the subtropical (∼ 30 to ∼ 40 • S), transitional (∼ 40 to ∼ 50 • S), subpolar (∼ 50 to ∼ 65 • S) and Antarctic (south of ∼ 65 • S) bands.Each of these biomes shows consistent ecological responses to 21st century climate change across most of the CMIP5 models studied.We further find that this banded structure is in general qualitative agreement with patterns and mechanisms of phytoplankton distribution changes which have emerged from observations over recent decades.

CMIP5 model description
A list of the models used along with relevant model details are summarized in Table 1.The scenarios used in our study are the historical and rcp8.5 scenarios from the IPCC's Fifth Assessment Report, with output data downloaded from http://cmip-pcmdi.llnl.gov/cmip5/data_portal.html.See Table 2 for a description of the variables downloaded and how they were used within this study.For all model analyses conducted here we use yearly time series, which were sometimes calculated from CMIP5 monthly output and sometimes taken straight from CMIP5 yearly output depending on availability.Some models lacked output data for certain variables.Table S1 shows which models had output for which variables.Only the first ensemble members (r1ip1) within the archives are used here.The historical scenario, spanning years 1850-2005, is forced with observed atmospheric CO 2 concentrations and is used to represent present-day conditions.The rcp8.5 scenario, spanning years 2006-2100, is representative of future unmitigated climate change conditions with radiative forcing increasing by 8.5 W m −2 relative to preindustrial by year 2100.See Taylor et al. (2012) and van Vuuren et al. (2011) for further details on CMIP5 experimental design and forcing scenarios.Absolute 100-year mean changes are calculated as the mean value from years 1980-1999 within the historical simulation subtracted from the mean value from years 2080-2099 within the rcp8.5 simulation.Relative change is defined as the 100-year absolute change divided by the historical 1980-1999 mean.
For multi-model statistical analysis, we weight models based on their similarity to avoid double counting and to preserve model independence.If two models are very similar in terms of their ocean biogeochemistry or physics (typically because they are two slightly different versions of the  same basic model coming from the same modelling centre -see Fig. S1 comparing phytoplankton biomass changes in HadGEM2-CC and HadGEM2-ES, for example), we give them each a weight of 0.5 instead of 1. See Table 1 for a list of model weights and Cabré et al. (2014) for a more detailed discussion on weighting.We do not attempt to weight models according to how well they reproduce observed chlorophyll a (chl) concentrations or primary productivities for the following reasons: (1) we cannot tell whether they reproduce current mean-state values of these variables for the right reasons, and (2) we would like to understand equally the reasons for each individual CMIP5 model's predictions and the reasons for the entire suite's predictions on average.

Bootstrap analysis (Figs. 1, 5, 6)
To quantify the significance of multi-model mean 100-year trends, we calculate the percentage of simulated model realizations that agree on the sign of a predicted trend for a given variable, using the statistical technique known as bootstrapping.We built 1000 realizations of the 100-year trend by randomly selecting n models (where n is the number of models with data available for any given variable) with replacement among the n available models.Within a single realization, one model may be represented more than once, while other models may not be represented at all.We take into account interannual variability by randomly selecting one of the 20 years from the present-day historical scenario (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) and one of the 20 years from the future rcp8.5 climate change scenario (2080-2099) for each selected model.For every variable of interest at every spatial grid point, we then create a realization of the 100-year trend by finding the difference between the two randomly chosen years.We then obtain the multi-model significance of this trend at each grid point by calculating the percentage of 1000 realizations that predict a positive change.Thus, the higher (lower) the bootstrap percentage above (below) 50 %, the greater the significance of the positive (negative) trend at a given location.This bootstrapping procedure provides a more robust measure of significance than simply calculating the percentage of models that agree based on single model runs alone because it both takes into account interannual variability and greatly increases the number of permutated realizations.See Cabré et al. (2014) for further details on application of the bootstrapping method to the CMIP5 data set.

Results and discussion
In this study, we attempt to understand how the general characteristics of SO phytoplankton may change with future warming by investigating biomass and productivity at both peak bloom times and averaged over the entire year.To this end, we choose to study the following two variables: (1) maximum annual surface phytoplankton biomass (henceforth PB, representative of phytoplankton biomass at the peak of an annual bloom) and (2) average annual primary production vertically integrated down to 100 m depth (henceforth PP, representative of average yearly water column integrated conditions).We conducted all of our analyses with both of these variables, but only show results for the variable which made the most sense to use in the context of the analysis.For example, whenever we analyse individual models, we show PB because we frequently only have monthly model output (with which to generate maximum, minimum or average annual data) at the surface of the ocean (i.e.monthly NO 3 , iron and light output are only available at the surface) and want to keep the variables we are cross-correlating spatially consistent whenever possible (either all variables at the surface only or all vertically integrated only).Although PB and PP are obviously different biological quantities (PB is surface phytoplankton biomass concentration and is directly affected by grazing, while PP is the integrated product of growth rate and biomass and is only indirectly affected by grazing -see references cited in Table 1 for model equation details), the direction of projected changes in the two variables are highly similar in our regions of interest (Fig. 1a,b;.Some exceptions to this occur between ∼ 50 and 65 • S in models GISS-E2-H-CC and CESM1-BGC (Figs.S1-S2); here PP increases while PB decreases, suggesting that the effects of top-down controls (grazing) win out over the effects of bottom-up controls (nutrients, light, temperature).Among the other models as well as other regions within these two models, however, changes in bottomup controls appear to explain most of the projected phytoplankton response such that patterns of predicted PP and PB change overlap significantly.Because of this and large uncertainties in how well the models' grazing parameterizations approximate the real ocean due to their incomplete food-web dynamics (see references cited in Table 1 for model equation details), we focus mostly on understanding the effects of bottom-up controls within all of the models.One other notable difference between PB and PP is that trends in PP appear to be slightly more regionally consistent across the models than trends in PB (Figs. , so that whenever we look at relationships across models, we use PP instead of PB.PP output is also available for a larger number of the models.

3.1
Zonally banded all-model mean 100-year changes (Fig. 1) Predicted multi-model mean 100-year changes in both PB and PP exhibit a zonally banded pattern similar to those predicted by individual models alone (Figs. 1a, b, S1-S2; Tables S2-S3 in Supplement).This leads to a natural division of the SO into four zonally banded biomes separated by switches in the sign of predicted PB and PP changes, as follows: 1. Subtropical -Within the first zonal band (∼ 30 to ∼ 40 • S), there is a predicted decrease in PB, PP, and wintertime nitrate concentrations (Figs.1c, S3).Here shallower wintertime MLDs (Fig. S4) and resulting decreases in nitrate supply are associated with increases in water column stratification and the climate-driven poleward expansion of subtropical gyres observed across all CMIP5 models (Meijers et al., 2012;Cabré et al., 2014).These abovementioned factors are proximate physical and biogeochemical drivers of predicted phytoplankton responses within the models, but what is the ultimate driver of all of these physical and biogeochemical changes?

Transitional
Historical and projected 21st century increases in the strength of the principal mode of variability in the SOcalled the Southern Annular Mode (SAM) -due to a combination of elevated CO 2 concentrations and ozone depletion could be one explanation.One highly agreed upon dynamical change captured within all of the CMIP5 models analysed here is an intensification and poleward shift of the SO westerly wind belt (Figs.1h, S10) associated with an increasingly positive phase of the SAM with future warming, as seen both here (Fig. S11) and in previous work (e.g.Yin, 2005;Arblaster and Meehl, 2006;Russell et al., 2006;Gillett and Fyfe, 2013;Zheng et al., 2013).This highly consistent increase in wind stress (which is most pronounced in the summer -plots not shown) south of 50 • S may explain the deepening of summertime MLDs south of 50 • S, while the decrease in wind stress between 30 and 50 • S may explain the shoaling of summertime MLDs in that region (Figs. 1d,S5).These changes in MLD can then affect nutrient supply to the surface, perhaps leading to the large decreases in surface nitrate concentrations between 30 and 50 • S (Figs. 1c, S3).
Warming, tropospheric stability changes, and southwardshifted storm tracks can also lead to shifts in cloudiness (e.g.Yin, 2005;Bender et al., 2012;Ceppi et al., 2014;Kay et al., 2014), which may help explain the increase in summertime cloud cover south of 50 • S (Figs. 1g, S8) and the concomitant decrease in summertime IPAR between 50 and 65 • S across the models (Figs.1f, S7).South of ∼ 65 • S in most models, sea ice melt (see Figs. S9; Turner et al., 2013;Mahlstein et al., 2013) allows more light to reach the surface of the ocean, resulting in a net increase in IPAR despite concurrent cloud cover increases.A robust analysis of the effects of SAM and SO westerly wind stress changes on the various proximate drivers we study here is out of the scope of this paper, but is a key issue that should be addressed in future work.
As for the ultimate driver of increases in surface iron concentrations, which contribute to increases in PB and PP in the Transitional (∼ 40-50 • S) and Antarctic (south of 65 • S) bands, there may be other complicating factors at work.Parameterizations of the marine iron cycle differ from model to model and include processes such as atmospheric dust deposition, phytoplankton-community dependent biological uptake and remineralization, vertical particle transport, scavenging, and the release of iron from sediments (e.g.J. K. Moore et al., 2013) is kept constant in the CMIP5 simulations, other processes listed above may change, thus altering surface iron inventories.For example, the increase in iron in the 40-50 • S Transitional band can be explained by enhanced vertical supply due to deeper wintertime mixed layers (Fig. S4) or by increases in summertime water column stratification, which can trap and concentrate iron deposited from the atmosphere closer to the surface.On the other hand, Misumi et al. (2014) showed that in the CESM1-BGC model (rcp8.5 scenario), a southward expansion of the subtropical gyre and changes in low-latitude iron utilization resulted in increased lateral advection of iron into the SO over the 21st century.Another potential ironenhancing mechanism in the SO is increased release of iron from sediments, a mechanism important within at least the GFDL models (J.Dunne, personal communication, 2014).

Interannual to decadal timescale analysis (Fig. 2)
To check the significance and robustness of the associations between phytoplankton abundance and the physicalbiogeochemical variables of interest (the bottom-up controls) discussed in Sect.3.1, we use regression and correlation analysis to study these associations in greater detail within three individual models with well-established, complex ocean biogeochemical modules (GFDL-ESM2G, HadGEM2-ES and IPSL-CM5A-MR).An important point to note is that multimodel mean changes in a given variable may be dominated by models with the biggest changes in some cases, so the analysis of individual models here is helpful in better illuminating particular relationships between variables.In Fig. 2, we show scatter plots of PB versus our variables of interest on multiple timescales, across the four chosen zonally banded biomes (as defined in Sect.3.1).We use only the grid points within each of the four zonal bands in a given model where the 100-year change in PB is predicted to go in the same direction as the entire band in the allmodel average.As an example, in Fig. 1a, we see that PB is expected to increase with future warming in the Antarctic band (south of 65 • S) in the all-model average; accordingly, we mask the grid points south of 65 • S within each individual model where PB increases and study those grid points alone to understand what is driving PB increases within the south-of-65 • S band as a whole.By this same procedure, we mask and investigate only the areas where PB decreases between 30-40 • S, where it increases between 40-50 • S, where it decreases between 50-65 • S, and where it increases south of 65 • S. We undertake this masking procedure because we would like to tease out the dominant driver of the net phytoplankton response within the zonal band of interest and masking helps to further amplify the signal we are looking for by focusing on what the majority of points we are interested in are doing, thus effectively diluting the con-founding effects of natural background variability.To confirm that masking does not significantly alter our results besides by potentially enhancing the signal-to-noise ratio of our correlations, we repeat any analyses (namely, Figs.2-4) involving masking with all (both masked and unmasked) grid points.Results from these all-inclusive analyses agree with those presented here for masked points only, but with slightly weaker correlation coefficients between phytoplankton biomass or productivity and a given driving variable of interest in some cases, as expected (see discussion below).
After spatially averaging PB and our variables of interest over the masked grid points within each zonal band, we then created different time series representing multiple timescales.To remove the effects of climate change and isolate interannual variability, we subtracted a 25-year running mean from every spatially averaged yearly data point within each scenario's raw yearly time series (historical from 1911-2005, rcp8.5 from 2006-2100).To capture variability and mechanisms which act on a longer than interannual but shorter than decadal timescale in the absence of confounding climate change effects, we took the 5-year running mean of the raw yearly time series data and then subtracted a 25year running mean from each 5-year running mean-smoothed annual value.Here we purposely chose to use detrended historical scenario time series rather than preindustrial control scenario time series (forced with constant preindustrial CO 2 concentrations) for practical reasons (not all the models provided all the necessary variables in the preindustrial control experiment).We did, however, prove that in at least model GFDL-ESM2G, the interannual drivers affect phytoplankton biomass in the same direction and with a similar magnitude in the preindustrial control case and the detrended historical and rcp8.5 cases, as expected.Finally, to investigate and emphasize the effects of climate change, we created a set of historical and rcp8.5 "climate change" time series by taking 10-year averages (not running means, but rather averages of non-overlapping 10-year intervals) of the same raw yearly spatially averaged time series as before.
A quick summary of the making of Fig. 2 is as follows: we spatially averaged PB and our driving variables of interest over each masked zonally banded biome and temporally correlated them across three distinct timescales of variability.Only those variables of interest which were significantly correlated with PB over at least two of the three (interannual, 5-year and decadal) studied timescales are shown (see Fig. S12 for examples of how correlations between PB and variables which were not chosen looked in comparison to the correlations between PB and the variable which was chosen).The driving variables shown are thus the ones whose relationships with PB hold on interannual, five-year, as well as longer-term climate change timescales in both the historical and rcp8.5 scenarios.Significantly, this implies that changes in these particular variables are the likely drivers of changes in phytoplankton biomass on an interannual as well as a longer-term timescale associated with future warming.2. * Wintertime MLD was also significant on all three timescales.* * Summertime MLD and avg annual cloud cover were also significant on all three timescales.
* * * Wintertime MLD was also significant on all three timescales.* * * * The y axis here is PP (µmol m −2 s −1 ) instead of PB because no variables were significantly correlated on at least two timescales with PB.Summertime IPAR was also significant on the same timescales as average annual sea ice cover.
It is important to note here that these inferred linkages are based only on correlations, but in all cases are also supported by model equations and previous studies.PB between 30 and 40 • S is strongly positively correlated with maximum annual surface nitrate concentration in all three models on all timescales (Fig. 2a).This suggests that predicted future decreases in PB between 30 and 40 • S are largely driven by climate warming-induced decreases in macronutrient supply to the surface during winter.This decreased supply is in turn a consequence of increased water column stratification and decreased maximum annual wintertime MLD associated with future warming, as was suggested by the analysis of multi-model mean maps discussed in Sect.3.1.
Between 40 and 50 • S, projections of enhanced PB are driven by either increases in iron concentrations (GFDL-ESM2G and IPSL-CM5A-MR) or reduced light limitation associated with shoaling of the summertime mixed layer (HadGEM2-ES) (Fig. 2b), again in agreement with the analysis in Sect.3.1.
Within the 50-65 • S band, where PB is predicted to decrease across all three models, light and iron are the most important limiting factors (Fig. 2c).For GFDL-ESM2G, in regions within this band where PB decreases with climate change, cloud cover (which is negatively correlated with PB -plot not shown) increases, leading to a concomitant decrease in surface light availability, which is positively correlated with PB (Fig. 2c).Furthermore, in both GFDL-ESM2G and HadGEM2-ES, the summertime MLD is predicted to deepen with climate change, creating an even more lightlimited environment for phytoplankton here (Fig. 2c), as was deduced using multi-model means in Sect. to the first two models and missing from the analysis in Sect.3.1, IPSL-CM5A-MR's wintertime surface iron concentrations appear to play the biggest role in determining PB between 50 and 65 • S (wintertime MLD was also well correlated with PB on all studied timescales here, most likely because it drives the supply of iron from the deep ocean) (Fig. 2c).South of 65 • S, iron is significantly positively correlated with PB on all three timescales within the models GFDL-ESM2G and IPSL-CM5A-MR (Fig. 2d), as was expected from multi-model mean change analyses in Sect.3.1.For HadGEM2-ES, both average annual sea ice fraction and maximum annual IPAR were significantly correlated with PP south of 65 • S, such that available light at the ocean surface is likely the limiting factor within this model's Antarctic band.An increase in IPAR due to a decrease in sea ice fraction is thus the most probable cause of projected phytoplankton abundance increases here, again in agreement with the reasoning in Sect.3.1.
In sum, these findings from Fig. 2 agree well with those deduced from Fig. 1, as discussed in Sect.3.1.In particular, within each zonally banded biome, the proposed drivers of projected phytoplankton responses in the all-model means are the same ones driving phytoplankton responses within the individual models studied here (with the possible exception of the 50-65 • S band, where iron appears to play a role within IPSL-CM5A-MR, but not in the all-model mean).Results for Fig. 2 were the same but with slightly smaller correlation coefficients when using all (both masked and unmasked) grid points (see Fig. S13).

Centennial timescale analysis (Fig. 3)
To confirm that the bottom-up controls on PB proposed in Sects.3.1 and 3.2.1 hold across the four SO biomes on even longer 100-year timescales, we undertake a spatial correlation analysis within the same three models as before.In Fig. 3, we show the results of this spatial correlation analysis in which we look at the relationship between 100-year changes in PB and the variables of interest at every grid point within each masked latitudinal band.Each dot in a scatter plot represents a masked grid point which undergoes a 100-year change in a variable of interest and an associated change in PB at that same grid point.By plotting only those variables with the largest magnitude correlation coefficients when correlated with PB, we are able to discover which variables affect PB most in each latitudinal band over 100-year timescales within each of the three models studied in detail (see Fig. S14 for an example of how correlations between PB and variables not chosen looked in comparison to correlations between PB and the variable chosen).For each chosen variable, scatter plots of either relative or absolute 100year changes are shown, depending on which type of change generated the clearest relationship between PB and the variable of interest.Least squares best-fit lines are drawn for each scatter plot to help visualize the slopes and enable comparison with the corresponding slopes in Fig. 2. Because it is difficult to accurately test for significance in this type of spatial correlation (neighbouring grid points are likely highly correlated, leading to large significance overestimates), these regression lines may or may not be statistically significant.Thus, the lines are meant only to serve as a qualitative visual guide.As in Fig. 2, Fig. 3 showed the same results but with potentially slightly smaller correlation coefficients when using all (both masked and unmasked) grid points (see Fig. S15).
Together, Figs. 2 and 3 show that the variables of interest within each zone that drive PB on decadal and shorter timescales also tend to be those that drive PB on an even longer 100-year climate change timescale.There are, however, a couple of important discrepancies.The first occurs in GFDL-ESM2G within the 50-65 • S band, where light limitation is shown to be most important on decadal and shorter timescales (Fig. 2c) while iron limitation appears to take over on a 100-year timescale (Fig. 3c).This suggests the presence of iron-light co-limitation in this region within GFDL-ESM2G, in agreement with previous studies (e.g.Sunda and Huntsman, 1997;Boyd et al., 2001;Feng et al., 2010).The second discrepancy occurs in IPSL-CM5A-MR within the south-of-65 • S band, where iron limitation is most important on decadal and shorter timescales (Fig. 2d) while increases in sea surface temperature (SST) become the dominant driver of PB increases on a 100-year timescale (Fig. 3d) (though iron is still somewhat important on the centennial timescale with R = 0.703 when spatially correlated with PB change -plot not shown).
In sum, we find that for the most part, the mechanisms within each zonal band that determine PB on decadal and shorter timescales tend to be those that determine PB on longer, centennial climate change-driven timescales as well.The magnitude of each driver's effect on phytoplankton biomass (as seen from the slopes of best-fit lines in Figs.2-3, summarized in Table 3) also remains the same across the relevant timescales, further supporting the notion that the same mechanisms act on the different timescales studied.Importantly, the magnitude of 100-year changes in the chosen variables of interest are also hypothetically large enough to drive most of the 100-year change in PB.We note, however, that in the real ocean, phytoplankton adaptation and evolution could alter the driver-response relationship observed at the interannual scale within these models.

Consistency of trends and mechanisms driving
phytoplankton changes across all models 3.3.1 Drivers of 100-year phytoplankton changes across all models (Fig. 4)  explicit phytoplankton biology.To answer this, we plot 100year changes in PP versus 100-year changes in chosen variables of interest across all of the models and look for amongmodel agreement as to the effect of these variables on PP within each masked zonal band (Fig. 4).These variables were chosen by first plotting all of the potential drivers of interest (listed in Table 2) and then choosing the ones which showed the strongest correlations or most consistent direc-  fraction and primary production south of 65 • S, we do not box the orange points in the PP versus average annual cloud cover plot (Fig. 4e) because we know that an increase in cloud fraction would decrease light availability and consequently lead to decreases, not increases, in primary production.Thus, we can safely ignore changes in cloud cover as a driver of changes in primary production among the models south of 65 • S and instead view these changes in cloud cover as merely a consequence of underlying dynamical changes already occurring in that region.Via this technique, we find a consistent set of mechanisms driving 100-year changes in productivity across the CMIP5 model suite (highlighted by coloured boxes in Fig. 4), in agreement with the mechanisms brought to light by the analyses in Sects.3.1 and 3.2.Nitrate emerges as the driver for changes in PP within the 30-40 • S band across all models (i.e.all red points lie in the third quadrant and within the red box in Fig. 4b).Models with greater relative decreases in wintertime surface nitrate concentrations undergo significantly (p < 0.05) greater decreases in average production within the 30-40 • S band.It is worth noting that this is the only significant, highly linear intermodel relationship within any of the zonal bands.In the rest of the bands, we mostly interpret only the sign, rather than the linearity, of PP changes related to the driving variables of interest across models.Within the 40-50 • S band, in general, models with increases in relative iron concentration and decreases in summertime MLD also experience relative increases in PP (Fig. 4a, c, purple boxes).Mod-els NorESM1-ME and IPSL-CM5A-LR are exceptions to this, however, in that PP still increases while iron concentrations decrease (Fig. 4c, purple unboxed).In these models, increases in light availability due to shoaling of summertime MLDs (Fig. 4a) and decreases in cloud cover (Fig. 4e) are large enough to cancel out the PP-suppressing effects of iron concentration decreases (Fig. 4c) between 40 and 50 • S. Further solidifying the importance of climate-driven changes in light availability within the 50-65 • S band, models predicting relative increases in summertime MLD or average annual cloud cover, along with decreases in maximum annual IPAR, also predict relative decreases in PP in this region (Fig. 4a,  d, e, green boxes).Iron also emerges as a potential driver of PP decreases within the 50-65 • S band, but not across all of the models (Fig. 4c).In models which undergo PP decreases concurrent with iron concentration increases (GISS-E2-R-CC, GISS-E2-H-CC, HadGEM2-CC, HadGEM2-ES, IPSL-CM5A-LR, NorESM1-ME, and MPI-ESM-LR; see Fig. 4c, green unboxed), reductions in light availability tend to be relatively large such that they win out in determining overall PP change.For example, GISS-E2-R-CC exhibits the largest relative iron increase between 50 and 65 • S out of all the models (Fig. 4c), but also the greatest relative summertime MLD deepening (Fig. 4a), leading to vast reductions in light supply to phytoplankton during the most productive time of year.An increase in both IPAR and iron supply across the models results in PP increases south of 65 • S, as highlighted by the orange boxes in Fig. 4c and d  LR, IPSL-CM5A-MR, and GFDL-ESM2G deviate from this trend slightly in that they experience small relative decreases in IPAR south of 65 • S, while still experiencing increases in PP.However, these three models also exhibit shoaling of the summertime MLD here, which would increase light availability, likely cancelling the effects of decreased IPAR at the surface.Note that for all models except for the three just mentioned, IPAR increases despite an increase in cloud cover (Fig. 4e, orange dots).This suggests that sea ice fraction, rather than cloud cover, is the most important factor in determining IPAR changes in this region within most models.
As sea ice cover declines near the Antarctic continent within the models (Fig. S9), more light is able to reach the ocean surface, ultimately leading to increased IPAR and PP here.
While general agreement on the mechanisms driving 100year phytoplankton changes among models is high, one noteworthy result is that there appear to be two distinct groups of models: one group with phytoplankton which are highly sensitive to changes in iron concentrations south of ∼ 40 • S (consisting of GFDL-ESM2, CESM1-BGC, IPSL-CM5A, CMCC-CESM -see Fig. S16, for models where zonal PB or PP changes closely follow zonal iron changes) and a second group with phytoplankton which are less iron sensitive (NorESM1-ME, HadGEM2, GISS-E2, MPI-ESM) or do not include iron at all (CanESM2, MIROC-ESM, MRI-ESM1).Iron cycling within the ocean remains poorly characterized and is typically crudely parameterized (if at all) compared to the macronutrients.These models also differ considerably in many aspects of their treatment of iron including but not limited to the magnitude and location of sources (from both the atmosphere and the sediments), ligand dynamics, scavenging losses, and iron to carbon biomass ratios (J.K. Moore et al., 2013).It is out of the scope of this paper to assess all of these differences, but at first glance, it appears that the models with more complex iron cycling dynamics have phytoplankton that are more sensitive to iron changes.For example, the more iron-sensitive GFDL-ESM2, CESM1-BGC, and IPSL-CM5A models have variable iron to carbon ratios and include sedimentary sources of iron (however crudely parame-terized) (Dunne et al., 2013;J. K. Moore et al., 2013;Aumont and Bopp, 2006), while the less iron-sensitive NorESM1-ME, HadGEM2, GISS-E2, MPI-ESM models do not (Assmann et al., 2010;Collins et al., 2011;Gregg, 2008).Models within the more iron-sensitive group tend to exhibit less well defined latitudinally banded 100-year phytoplankton changes, while the other models tend to exhibit a more obviously banded PB and PP change structure (see Figs. S1-2).These less iron-sensitive models also frequently display iron and phytoplankton changes of opposite signs south of 40 • S (Fig. 4c, unboxed purple and green points; Fig. S16).In these cases, changes in light availability due to changes in MLD and IPAR are able to explain predicted phytoplankton trends (see Fig. S16, for models where zonal PB or PP changes closely follow zonal MLD and/or IPAR changes).Within the group of models with iron-sensitive phytoplankton, changes in physical variables altering light availability are also occurring, but their effects are much less pronounced because iron plays a more dominant role.As was discussed in Sect.3.1, changes in MLD and IPAR in both groups of models are in turn driven by first-order changes in oceanatmosphere dynamics associated with climate warming and an increasingly positive SAM index, such as westerly wind intensification, alterations to tropospheric stability and thermal structure (e.g.Ceppi et al., 2014;Kay et al., 2014), and poleward displacement of extratropical storm-tracks and associated clouds (e.g.Yin, 2005;Bender et al., 2012).

Spatial agreement on projected changes across all models (Figs. 5-6)
To get a wider sense of spatial agreement among models throughout the SO, we look at maps of intermodel consistency in projected SO phytoplankton trends and their proposed drivers across all 16 CMIP5 models with ocean biogeochemistry in Fig. 5 (complementary to Fig. 1).The maps in Fig. 5 detail the fraction of model realizations (via the bootstrap technique explained in Sect.2.2) that predict a positive trend in the listed variable at each grid point.Thus, the closer the fraction to 1 at a given location, the greater the intermodel agreement on a positive trend at that point, and the closer the fraction to 0, the larger the intermodel agreement on a negative trend at that point (0.5 denotes the greatest amount of intermodel disagreement, where 50 % of model realizations predict an increasing trend and 50 % predict a decreasing trend).To also get a better idea of how well models agree with one another within each zonal band, Fig. 6 shows zonally averaged all-model mean projected trends (zonal averages of Fig. 1) and zonal band averaged intermodel agreement percentages (areal averages over each zonal band of Fig. 5, listed above each zonal band accompanied by an arrow indicating the direction of the trend agreed upon by the majority of model realizations).Only percentages for variables which are most important within each zonal band (as determined by Figs.1-4) are listed and as such, represent a summary of the important drivers of projected phytoplankton change discussed here.Agreement among models is highest at the centre of each zonal band (Fig. 5), but decreases towards the edges due to offsets in the precise boundaries of water masses among the models.These slight offsets lower the zonal band-average agreement among models shown in Fig. 6, such that if one were able to perfectly compare water masses among models, consistency in predicted trends within each zonally banded biome would likely be even higher.Within the subtropical (30 to 40 • S) band, the majority of model realizations predict a decrease in both PB (64 %) and PP (62 %), accompanied by a highly consistent decrease in wintertime nitrate supply (77 %) (Figs.5a-c, 6).This projected change agrees with the general expectation from previous theoretical and modelling studies that warming should stratify the water column, decrease macronutrient supply, and consequently lower biological productivity within the subtropics (e.g.Sarmiento et al., 2004b;Doney, 2006;Cabré et al., 2014).Within the transitional (40 • S to 50 • S) band, most of the model realizations predict an increase in PP (70 %) while only around half of the models predict an increase in PB (55 %) (Figs.5a-b, 6).These predicted shifts are accompanied by a decrease in summertime MLD (71 %) and an increase in wintertime iron concentration (64 %) (Fig. 5cd; Fig. 6).Because of a predicted poleward shift of the westerly winds in all of the models (Figs.5h, S10), winds will weaken here, shoaling the MLD and prolonging the growing season by allowing phytoplankton to remain within the well-lit surface layers for longer.Thus, enhanced future phytoplankton populations within this transitional band are not unexpected.Within the subpolar (50 to 65 • S) band, mod-Figure 6. Summary of predicted phytoplankton responses and their drivers within each zonal band.Here each coloured line represents normalized 100-year all-model mean zonal changes in the listed variable.Each variable was normalized by first computing the all-model mean zonally averaged 100-year change at every latitude and then dividing by the absolute value of the largest of these changes occurring south of 30 • S. Listed above each band is the spatially averaged percentage of model realizations that agree on the sign of the change (based on a bootstrap analysis -see Sect.2.2 and Fig. 5) in each variable over that band.The coloured arrows denote the direction of the trend agreed upon by the majority of models.The number of models (n) and the total model weight (w) taken into account for each variable are listed in Fig. 5. els are not as consistent in their predictions of phytoplankton changes compared with the other regions.The majority of model realizations predict an increase in PP (59 %), while 55 % predict a decrease in PB (Fig. 5a-b).Predicted changes in driving variables are somewhat more consistent within this region, however, with a decrease in summertime MLD predicted by 56 % of model realizations, a decrease in IPAR predicted by 71 %, and an increase in cloud cover predicted by 60 % (Figs.5c, f, g, 6).With a projected poleward shift of the westerlies, cloud cover should increase (due to a concomitant shift in storm-track cloudiness and/or altered tropospheric stability with future warming) and MLDs should deepen as winds intensify within this band, both of which act to decrease phytoplankton populations, exactly as we see here.Within the Antarctic (south of 65 • S) band, 76 % of model realizations predict an increase in PP, while 64 % predict an increase in PB, both of which are associated with projected increases in wintertime iron concentrations (72 %) and summertime light availability (59 %) 6).This goes with our expectation that the melting of sea ice projected by the models will lead to higher amounts of light reaching the ocean surface and that intensified westerlies will bring a larger supply of upwelled iron to the surface in this region, both of which act to increase phytoplankton populations, just as we see here.

Linking CMIP5 model projections to observations
Because the same interannual mechanisms for phytoplankton growth hold on 5-year, decadal, and even longer-term timescales within the CMIP5 models, it is reasonable to compare recent observations to future model projections if it is also assumed that short-term drivers of observed phytoplankton variability propagate up to longer-term timescales in the real ocean as well.However, it is out of the scope of this paper to compare recent observations to historical model output from the same period.Instead, we would like to understand how our modelled 21st century SO predictions compare to observed mechanisms and trends thus far.
The SO satellite chl record is not yet long enough to separate the effects of climate change from those of interannual processes driven by the leading modes of shorter-term variability in the SO (e.g.Boyd et al., 2008;Strutton et al., 2012;Henson et al., 2010;Beaulieu et al., 2013), the most important being the SAM (Thompson and Solomon, 2002); in situ data from field campaigns suffer from the same tem-poral constraint.Furthermore, while models generate perfectly continuous data, observations tend to contain many more gaps, such that a longer observational time series is needed to detect significant trends compared to model data when the same threshold of significance is applied.For these reasons, many observational studies have looked at the effects of SAM and other modes of variability, rather than climate change, on phytoplankton abundance and productivity.These types of studies can, however, still provide essential insight into the mechanisms driving possible longer-term changes.For example, as was mentioned before, the SAM index is expected to become increasingly positive as SO westerlies strengthen and move poleward with future warming (see Fig. S11 for projected SAM time series within the CMIP5 models).We have shown here that at least within the CMIP5 models, mechanisms responsible for changes in phytoplankton biomass on interannual and five-year timescales are also responsible for projected 100-year trends within the SO.Thus, understanding the effects of a more positive SAM on SO phytoplankton may help predict the direction of phytoplankton changes in a warmer future climate.
Another important caveat to keep in mind when looking at observational data is that observations rarely span consistent time frames, making it difficult to compare studies in a perfectly congruent way.For instance, it has been shown that the magnitude and sign of observed trends can be very sensitive to the start and end years analysed (e.g.Fay et al., 2014).Thus, rather than directly comparing recently observed trends with 21st century CMIP5 projections, we seek only to qualitatively understand whether there are common mechanisms and directions of change within the observational data and model projections.The observational studies cited in the following paragraphs are visually and tabularly summarized in Fig. 7 and Table S4.
Analysing satellite data over years 1997-2004, Lovenduski and Gruber (2005) (LG2005) found a negative correlation (though not significant) between SAM and chl concentrations within the SO Subtropical Zone (∼ 30-40 • S), due to increased stratification and decreased upwelling of macronutrients during positive SAM periods.Assuming that SAM will continue to increase with future warming and that the same driving mechanisms will hold on timescales ranging from interannual to centennial, phytoplankton biomass would be expected to decrease over the 21st century within the Subtropical Zone (∼ 30-40 • S) due to enhanced macronutrient limitation, which is indeed what the CMIP5 models predict.
Via a combination of satellite, reanalysis and model data, Johnston and Gabric (2011) (JG2011) found that both summertime chl concentrations and primary productivity increased within the Australian sector between 40 and 50 • S over the years 1997-2007, which they attribute to increased water column stratification or enhanced mineral dust deposition from Australia.Gregg et al. (2005) (G2005) likewise found an increase in chl concentrations just south of Australia (∼ 35-55 • S) from satellite data over the period 1998-2003, accompanied by an increase in springtime SST, likely associated with a shoaling of the mixed layer.Using satellite chl concentrations (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) calculated in two different ways, Siegel et al. (2013) (S2013) also reported an increase in chl concentrations between ∼ 40 and 50 • S.These proposed mechanisms and directions of trends are consistent with those of the CMIP5 models, which predict that increased dissolved iron concentrations together with decreased light limitation due to shallower MLDs during blooms will drive 21st century phytoplankton increases within the 40-50 • S band.
From satellite data (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004), LG2005 found a significant negative correlation between SAM and chl concentrations within the ∼ mid-40 to mid-50 • S latitudes below Australia, which they ascribe to increased light limitation due to deeper summertime mixed layers in positive SAM phases.Consistent with LG2005 and an increasingly positive SAM index, Takao et al. (2012) (T2012) found a decreasing trend in summertime net primary productivity within the Indian Ocean sector of the Polar Frontal Zone (centred slightly north of ∼ 55 • S) using satellite ocean colour data from 1997-2007.Within the Australian sector, JG2011 observed similar decreases in summer and springtime chl concentrations between 55-60 • S from 1997-2007, allegedly due to a decrease in the northward Ekman transport and supply of iron here.Based on both in situ shipboard measurements and satellite-derived chl concentrations, Montes-Hugo et al. (2009) (MH2009) also reported a decrease in phytoplankton biomass between 1978-1986 and 1998-2006 within the northern subregion of the West Antarctic Peninsula (61.8 to 64.5 • S) because of deeper summertime mixed layers, in turn driven by stronger winds and decreased sea ice extent.Compiling net haul data from nine different countries, Atkinson et al. (2004) (A2004) found significant decreases in krill density between 1976 and 2003 within the southwest Atlantic sector of the SO between ∼ 50 and 65 • S, which they attributed to decreases in phytoplankton populations.These findings fit with the previously discussed CMIP5 model predictions of 21st century decreases in phytoplankton biomass between ∼ 50 and 65 • S, which we attribute to more stressful light (as in LG2005 and MH2009) and/or iron conditions for phytoplankton (as in JG2011).
In the Antarctic Zone (south of ∼ 60 • S), Ayers and Strutton (2013) (AS2013) found a correlation between a more positive SAM and increased upwelling of nutrients based on multiple repeat hydrographic sections.LG2005 found a similar positive correlation between SAM and chl concentrations here due to increased upwelling and iron supply in positive SAM periods (also in agreement with a modelling study by Hauck et al., 2013).Again, assuming that SAM will continue to increase with future warming and that the same driving mechanisms will hold on timescales ranging from interannual to centennial, we expect increases in iron supply to drive phytoplankton biomass increases south of ∼ 60 • S Hatching denotes regions where trends calculated as least-squares best-fit lines to the time series are significant using a two-tailed t test at p < 0.05.with future warming, which is indeed what the CMIP5 models predict.In terms of trends, MH2009 report an increase in southern West Antarctic Peninsula (63.8 to 67.8 • S) summertime phytoplankton populations between 1978-1986 and 1998-2006, which they ascribe to decreased light limitation, driven by a decrease in cloudiness and wind intensity and an increase in the number of ice-free summer days.Meanwhile, S2013 observed a thin band of chl increase around ∼ 65 • S over the years 1997-2010.These observations are also consistent with future CMIP5 model projections, which predict that decreased sea ice cover will drive phytoplankton abundance increases south of ∼ 65 • S in spite of an increase in cloud cover (contrary to the decrease in cloudiness measured by MH2009).Lastly, Smith and Comiso (2008) (SC2008) calculate an increase in annual primary productivity over the entire Southern Ocean (defined as south of 60 • S) between 1997and 2006, while Arrigo et al. (2008) ) (A2008) calculate no significant trend over the same period.The discrepancy between these two works is partly due to the fact that A2008 define the Southern Ocean as south of 50 • S instead of 60 • S, and the region in between 50-60 • S underwent a decrease in productivity over both studies' time periods (reducing the magnitude of the increasing trend over the rest of the SO), again consistent with model projections of fu-ture phytoplankton biomass decrease between 50-65 • S and increase south of 65 • S.
To conduct our own analysis, we obtained monthly global satellite chl fields generated by the latest version of Sea-WiFS' (Sea-viewing Wide Field-of-view Sensor) band-ratio algorithm (OC4v6) (http://oceancolor.gsfc.nasa.gov/cgi/l3)from September 1997 to December 2010.The linear trend in Fig. 7c was calculated from yearly-averaged monthly chl anomalies, which ensures minimal autocorrelation.To look at trends in observed summertime MLD, monthly ocean temperature and salinity reanalysis products from the Met Office Hadley Centre's EN3 data set (http://www.metoffice.gov.uk/hadobs/en3/) were used to calculate minimum monthly MLDs for each year from 1950 to 2013.To look at trends in observed summertime cloud cover, synoptic monthly mean ERA-INTERIM (http://www.ecmwf.int/en/forecasts/datasets/era-interim-dataset-january-1979-present) reanalysis products of total cloud cover from December 1980 to February 2013 were averaged over the summer months (December to February) of each year to generate a yearly summertime cloud cover time series.
We found that recently observed spatial distributions of SO chl trends over the SeaWiFS period (1998-2010) (Fig. 7b, c In sum, the observed spatial distribution in trends of phytoplankton productivity, MLD and cloud cover over the past few decades qualitatively matches the latitudinally banded structure of the respective 100-year 21st century trends predicted by the CMIP5 models.We have found that (a) in CMIP5 simulations, interannual effects propagate up to 100-year timescales and (b) drivers for short-term biomass change are similar in models and observations within individual zonally banded biomes.If the CMIP5 model mechanisms and projections are to be trusted, then this suggests that observations may already contain a climate change signal even though this signal cannot be teased apart from decadal variability and shorter-term noise just yet (e.g.Henson et al., 2010).In agreement with discussions above, the fact that long-term model projections appear to agree with the sign of observed SAM-driven effects throughout the SO further suggests that an increasingly positive SAM may be responsible for the predicted zonally banded pattern of phytoplankton biomass changes in the models, though further work is needed to precisely quantify SAM's contribution to PB and PP changes and variability within the CMIP5 model suite.

Conclusions
The 16 CMIP5 models with explicit phytoplankton ecology predict a zonally banded pattern of 21st century phytoplankton biomass and productivity changes within the Southern Ocean: a decrease in the subtropical band (∼ 30-40 • S), an increase in the transitional band (∼ 40-50 • S), a decrease in the subpolar band (∼ 50-65 • S), and an increase in the Antarctic band (south of ∼ 65 • S).In line with previous studies, light (controlled by cloud cover, summertime MLD during blooms, and sea ice fraction) and iron supply are found to be the most important factors driving phytoplankton changes in the transitional and subpolar Southern Ocean (south of ∼ 40 • S), while nitrate is found to be the most important driving factor in the subtropical Southern Ocean (∼ 30-40 • S).Shifts in these driving variables consistently bring about changes in phytoplankton abundance and production on multiple timescales.In particular, within a given zonally banded biome in an individual model, the same mechanisms are generally responsible for phytoplankton biomass changes on an interannual, decadal and 100-year basis.This suggests that the mechanisms affecting shorter-term phytoplankton vari-ability, which can in principle be gauged from in situ or satellite observations, are also likely to be the mechanisms responsible for climate-driven phytoplankton changes over the 21st century.It is important to note that the relationships between phytoplankton responses and their potential drivers discussed here are based on correlative analysis and thus do not perfectly prove causation.It is promising, however, that in all cases the significant and most strongly correlated phytoplankton and potential driver relationships matched with expectations based on both previous studies and model equations.
Twenty-first century trends in phytoplankton productivity predicted by the CMIP5 models go in the same direction as observed trends over the last couple of decades and tentatively agree with the sign of established SAM-driven changes.This suggests that an increasingly positive SAM may be responsible for the projected zonally banded trends in phytoplankton productivity and biomass that we observe in the CMIP5 models, though more work is needed to carefully test this hypothesis.Additionally, since the observed trends in, and drivers of, short-term biomass change seem to agree with those of modelled decadal and centennial projections, it is possible that climate change is already having an effect on SO phytoplankton biology within the real ocean.
With such short and discontinuous observational records, our model-observational data intercomparison is clearly only qualitative at this point in time.We advocate for longer and more continuous in situ phytoplankton biomass and satellite chl data collection in this important but massively undersampled region of the ocean to allow for the emergence of a climate change signal from short-term variability.The main result of this study -a consistency of the model-projected phytoplankton trends within four distinct SO bands over the 21st century -suggests a framework for selecting a minimum number of sites for future SO biogeochemical observational time series stations or repeat sampling campaigns; at a minimum, one or two representative time series are needed from each of the four SO bands described here.These data sets (and any observational data sets, for that matter) are subject to all manner of spatial and temporal caveats, but over time and in combination with larger-scale satellite observations, longer-term in situ time series will allow us to distinguish natural variability from the climate change signal and more readily compare observed mechanisms and trends with those predicted by our models.
Follow-up work is needed to determine how projected changes in phytoplankton biomass and productivity will affect SO carbon and nutrient cycling, as well as how changes in the characteristics of regional SO seasonality can affect these long-term trends (Thomalla et al., 2011).Driving higher trophic-level models with projected CMIP5 phytoplankton abundances may also yield important insights into how ecologically and economically important species such as zooplankton, krill, marine mammals, penguins, and seabirds will respond to climate change.Given the critical importance of the SO in driving global carbon and nutrient cycles as well as low-latitude productivity, our results highlight the need for both long-term in situ and satellite monitoring of Southern Ocean biology and biogeochemistry.
The Supplement related to this article is available online at doi:10.5194/bg-12-5715-2015-supplement.

Figure 1 .
Figure 1.All-model mean 100-year changes.100-year all-model mean changes in (a) maximum annual surface phytoplankton biomass (PB), (b) average annual 100 m depth vertically integrated primary production (PP), (c) wintertime surface nitrate concentration, (d) summertime mixed layer depth (MLD), (e) wintertime surface dissolved iron concentration, (f) summertime incident photosynthetically available radiation (IPAR), (g) summertime percentage area of grid cell covered by clouds, and (h) average annual zonal wind stress.Hatched areas are where greater than 80 % of model realizations agree on the sign of the change using a bootstrap significance test (see Sect. 2.2 for methodological details).Zero contours for PP change are plotted over each map.The number of models (n) and the total model weight (w) taken into account for each variable are listed in Fig. 5. Historical all-model mean maps are presented in Figs.S1-10.

Figure 2 .
Figure 2. Drivers of phytoplankton biomass on multiple timescales.Scatter plots of PB versus the listed variable on interannual and 5year (both with their climate change signals removed) as well as 10-year timescales.Each column corresponds to a different model, while each row corresponds to a different zonal band.Slopes of the historical and rcp8.5 10-year average best-fit lines are listed.Only variables with significant (p < 0.05) best-fit lines on at least two out of the three timescales studied (interannual, 5-year, and 10-year) are shown.Best-fit lines are drawn only when correlations are significant (p < 0.05).Variables tested are all those listed in Table2.* Wintertime MLD was also significant on all three timescales.* * Summertime MLD and avg annual cloud cover were also significant on all three timescales.** * Wintertime MLD was also significant on all three timescales.* * * * The y axis here is PP (µmol m −2 s −1 ) instead of PB because no variables were significantly correlated on at least two timescales with PB.Summertime IPAR was also significant on the same timescales as average annual sea ice cover.

Figure 3 .
Figure 3. Spatial correlation scatter plots of 100-year changes in phytoplankton biomass versus 100-year changes in driving variables of interest.Each column corresponds to a different model, while each row corresponds to a different zonal band.Only the variable of interest with the largest magnitude correlation coefficient is plotted for each zone within each model.Variables tested are all those listed in Table2.Relative changes are plotted for PB vs. nitrate, while absolute changes are plotted for PB vs. all other variables.Best-fit lines are forced to have a zero-intercept.Correlation coefficients and slopes of best-fit lines corresponding to absolute 100-year changes (even for nitrate) to facilitate comparison with slopes in Fig.2are listed beneath the variable names.
tions of changes across the models, guided by the relationships found in Figs.1-3.Here each point in a scatter plot represents a 100-year change in PP versus a 100-year change in the variable of interest spatially averaged over the given model's masked zonal band.We box only the points driven by processes which could be logically predicted from previously discussed mechanisms or model equations.For example, although almost all models undergo increases in cloud

Figure 4 .
Figure 4. Drivers of 100-year phytoplankton productivity changes across CMIP5 models, by model and latitudinal band.Scatter plots of each model's 100-year relative change in PP versus its corresponding relative change in the listed variable within each zonal band.Each colour represents a different zonal band, while each symbol represents a different model.Coloured boxes enclose points which behave in line with our expectations and proposed mechanisms based on Figs.1-3.Best-fit lines are drawn only when correlations are significant (p < 0.05).

Figure 5 .
Figure5.Spatial agreement among models on the sign of predicted trends, represented by maps of the fraction of model realizations that agree on a positive 100-year change in variables of interest at each grid point, based on a bootstrap analysis test (see Sect. 2.2).The closer to 1 the grid point, the greater the agreement among models on an increase.The closer to 0 the grid point, the greater the agreement among models on a decrease.
Fig-ure S17 complements Fig. 6 by showing zonally averaged all-model historical means and 100-year absolute changes in all variables of interest.

Figure 7 .
Figure 7. Observed phytoplankton trends and variability.(a) Summary of past studies looking at trends and SAM-driven variability in phytoplankton biomass and productivity.Orange/red regions are areas where past studies have found positive trends in phytoplankton biomass or productivity, whereas blue regions are areas where past studies have found negative trends.Each coloured region or point is labelled with the corresponding publication.See Table S4 for further details on each study.(b) Average monthly SeaWiFS chl concentrations, along with yearly trends in (c) chl, (d) summertime cloud cover from ERA-INTERIM reanalysis, and (e) summertime MLD from Hadley reanalysis.Hatching denotes regions where trends calculated as least-squares best-fit lines to the time series are significant using a two-tailed t test at p < 0.05.
) generally correspond well with CMIP5 all-model mean pro-.1a, b), with the largest observed chl increases occurring between ∼ 40 and 50 • S and south of ∼ 65 • S and decreases occurring between ∼ 50 and 65 • S. We also found that spatial distributions of recent trends in summertime MLD and cloud cover generally match with CMIP5 model projections as well.For example, the largest observed increases in summertime MLD (over the years 1950-2013) and cloud cover (over the years 1980-2013) occur south of ∼ 50 • S, while the largest decreases occur north of ∼ 50 • S (Fig.7d, e compared with Fig.1d, g, respectively).

Leung et al.: A latitudinally banded phytoplankton response 5717Table 1 .
CMIP5 model details.Summary of all the CMIP5 models that keep track of phytoplankton biomass and/or primary production with information on the following for each model: spatial resolution in the atmosphere and ocean, explicitly modelled nutrients, ecology subroutine, references and weight (wt) applied in the all-model averages.( * Note: CMCC-CESM runs did not appear to reach equilibrium, which was a necessary condition we imposed in order to work with a model; thus, we only show data from CMCC-CESM in the Supplement figures, but do not take this model into account in the all-model averages.) Biogeosciences, 12, 5715-5734, 2015 www.biogeosciences.net/12/5715/2015/S.

Table 2 .
Details of the CMIP5 variables studied here.

Table 3 .
Summary of the best-fit line slopes between phytoplankton biomass and given variables of interest, corresponding to Figs.2 and 3.The first of the three numbers in each table entry is the historical 10-year average slope between PB and the given variable of interest, while the second number is the rcp8.5 10-year average slope, both from Fig.2.The third number in each entry is the 100-year change spatial correlation slope from Fig.3.For variable units, see Figs.2 and 3.