Articles | Volume 17, issue 18
https://doi.org/10.5194/bg-17-4633-2020
https://doi.org/10.5194/bg-17-4633-2020
Research article
 | Highlight paper
 | 
25 Sep 2020
Research article | Highlight paper |  | 25 Sep 2020

Increase in ocean acidity variability and extremes under increasing atmospheric CO2

Friedrich A. Burger, Jasmin G. John, and Thomas L. Frölicher
Abstract

Ocean acidity extreme events are short-term periods of relatively high [H+] concentrations. The uptake of anthropogenic CO2 emissions by the ocean is expected to lead to more frequent and intense ocean acidity extreme events, not only due to changes in the long-term mean but also due to changes in short-term variability. Here, we use daily mean output from a five-member ensemble simulation of a comprehensive Earth system model under low- and high-CO2-emission scenarios to quantify historical and future changes in ocean acidity extreme events. When defining extremes relative to a fixed preindustrial baseline, the projected increase in mean [H+] causes the entire surface ocean to reach a near-permanent acidity extreme state by 2030 under both the low- and high-CO2-emission scenarios. When defining extremes relative to a shifting baseline (i.e., neglecting the changes in mean [H+]), ocean acidity extremes are also projected to increase because of the simulated increase in [H+] variability; e.g., the number of days with extremely high surface [H+] conditions is projected to increase by a factor of 14 by the end of the 21st century under the high-CO2-emission scenario relative to preindustrial levels. Furthermore, the duration of individual extreme events is projected to triple, and the maximal intensity and the volume extent in the upper 200 m are projected to quintuple. Similar changes are projected in the thermocline. Under the low-emission scenario, the increases in ocean acidity extreme-event characteristics are substantially reduced. At the surface, the increases in [H+] variability are mainly driven by increases in [H+] seasonality, whereas changes in thermocline [H+] variability are more influenced by interannual variability. Increases in [H+] variability arise predominantly from increases in the sensitivity of [H+] to variations in its drivers (i.e., carbon, alkalinity, and temperature) due to the increase in oceanic anthropogenic carbon. The projected increase in [H+] variability and extremes may enhance the risk of detrimental impacts on marine organisms, especially for those that are adapted to a more stable environment.

Dates
1 Introduction

Since the beginning of the industrial revolution, the ocean has absorbed about a quarter of the carbon dioxide (CO2) released by human activities through burning fossil fuel and altering land use (Friedlingstein et al.2019). Oceanic uptake of anthropogenic CO2 slows global warming by reducing atmospheric CO2 but also leads to major changes in the chemical composition of seawater through acidification (Gattuso and Buddemeier2000; Caldeira and Wickett2003; Orr et al.2005; Doney et al.2009). When CO2 dissolves in seawater, it forms carbonic acid that dissociates into bicarbonate ([HCO3-]), releasing hydrogen ions ([H+]) and thereby reducing pH (pH =-log([H+])). The rise in [H+] is partially buffered by the conversion of carbonate ions ([CO32-]) to [HCO3-]. The associated decline in [CO32-] reduces the calcium carbonate saturation state Ω=[Ca2+][CO32-]/[Ca2+][CO32-]sat, i.e., the product of calcium and carbonate ion concentrations relative to the product at saturation. Undersaturated waters with Ω<1 are corrosive for calcium carbonate minerals. Each type of calcium carbonate mineral has its individual saturation state Ω due to different solubilities, e.g, ΩC for calcite and ΩA for aragonite. Over the last four decades the surface ocean pH has declined by about 0.02 pH units per decade (Bindoff et al.2020). Continued CO2 uptake by the ocean will further exacerbate ocean acidification in the near future (Caldeira and Wickett2003; Orr et al.2005; Bindoff et al.2020; Terhaar et al.2020), with potentially major consequences for marine life (Doney et al.2009) and ocean biogeochemical cycling (Gehlen et al.2012).

Superimposed onto the long-term decadal- to centennial-scale ocean acidification trend are short-term extreme events on daily to monthly timescales during which ocean pH and Ω are much lower than usual (Hofmann et al.2011; Joint et al.2011; Hauri et al.2013). These events can be driven by different processes, such as ocean mixing, biological production and remineralization, mineral dissolution, temperature and air–sea gas exchange variations, or a combination thereof (Lauvset et al.2020). In eastern boundary upwelling systems, for example, short-term upwelling events and mesoscale processes can lead to low-surface-pH events and to short-term shoaling of the saturation horizon (i.e., the depth between the supersaturated upper ocean and the undersaturated deep ocean; Feely et al.2008; Leinweber and Gruber2013). Ocean pH can also rapidly change as a consequence of microbial activity (Joint et al.2011). Phytoplankton blooms and accompanying respiration drastically increase the partial pressure of CO2 (pCO2) and reduce pH in the thermocline (Sarmiento and Gruber2006). Such extreme events may have pH levels that are much lower than the mean pH conditions projected for the near future (Hofmann et al.2011).

Most of the scientific literature on ocean acidification has focused on gradual changes in the mean state in ocean chemistry (Orr et al.2005; Bopp et al.2013; Frölicher et al.2016; Terhaar et al.2019b). However, to understand the full consequences of ocean acidification on marine organisms and ecosystem services, it is also necessary to understand how variability and extremes in ocean acidity change under increasing atmospheric CO2 (Kroeker et al.2020). The ability of marine organisms and ecosystems to adapt to ocean acidification may depend on whether the species have evolved in a chemically stable or a highly variable environment (Rivest et al.2017; Cornwall et al.2020). Furthermore, if the frequency and intensity of short-term extreme events strongly increase, in addition to the long-term acidification, some organisms may have difficulties in adapting, especially if key CO2 system variables cross some critical thresholds, e.g., from calcium carbonate supersaturation to undersaturation. Key plankton species such as coccolithophores (Riebesell et al.2000), foraminifera, and pteropods (Bednaršek et al.2012) were found to be adversely affected by low carbonate ion concentrations. After only several days of being exposed to waters which are undersaturated with respect to aragonite, some species such as pteropods already show reduced calcification, growth, and survival rates (Kroeker et al.2013; Bednaršek et al.2014). Carbonate system variability also plays a role in shaping the diversity and biomass of benthic communities (Hall-Spencer et al.2008; Kroeker et al.2011). In laboratory experiments in which deep-water corals are exposed to low-pH waters for a week, some corals exhibit reduced calcification, while recovery may be possible when the low-pH condition persists for several months, stressing the importance of high-frequency variability and short-term acidification events (Form and Riebesell2012). There is also growing evidence that the organism response to variability in ocean acidity could change with ocean acidification (Britton et al.2016). Therefore, understanding the temporal variability of ocean carbonate chemistry and how that will change is important for understanding the impacts of ocean acidification on marine organisms and ecosystems (Hofmann et al.2011).

Changes in extremes can arise from changes in the mean, variability, or shape of the probability distribution (Coles2001). There exists no general accepted definition of an extreme event beyond the common understanding that an extreme is rare (Weyer2019). As a result, many different approaches exist to define extreme events (Smith2011). If a relative threshold (e.g., quantile) is used to define an extreme event, it is important to distinguish between extreme events that are defined with respect to a fixed reference period or baseline, or if the reference period or baseline moves with time. If the baseline is fixed, the changes in the mean state as well as changes in variability and higher moments of the distribution contribute to changes in extreme events (e.g., Fischer and Knutti2015; Frölicher et al.2018; Oliver et al.2018). However, if a shifting baseline is used, changes in the mean state do not contribute to changes in extreme events (e.g., Stephenson2008; Seneviratne et al.2012; Zscheischler and Seneviratne2017; Cheung and Frölicher2020; Vogel et al.2020). In this case, changes in extremes arise solely due to changes in variability and higher moments of the distribution (Oliver et al.2019). This latter definition ensures that values are not considered extreme solely because the baseline changes under climate change (Jacox2019; Oliver et al.2019). Whether extreme events should be defined with respect to a fixed baseline or with respect to a shifting baseline depends on the scientific question. For example, the shifting-baseline approach may be more appropriate when the ecosystems under consideration are likely able to adapt to the mean changes but not to changes in variability (Seneviratne et al.2012; Oliver et al.2019). Here, we use both approaches, with a special focus on the analysis of ocean acidity extremes with respect to shifting baselines.

Under continued long-term ocean acidification (i.e., changes in the mean), one can expect that extreme events in [H+] and Ω, when defined with respect to a fixed reference period or baseline, will become more frequent and intense (Hauri et al.2013). In addition to the changes in the mean, recent studies suggest that the seasonal cycles in [H+] and Ω are also strongly modulated under elevated atmospheric CO2. Higher background concentrations of dissolved inorganic carbon and warmer temperatures produce stronger departures from mean state values for a given change in pertinent physical or chemical drivers for [H+] and weaker departures for Ω (Kwiatkowski and Orr2018; Fassbender et al.2018). Other studies have also addressed the changes in the seasonal cycle of pCO2 (Landschützer et al.2018; Gallego et al.2018; McNeil and Sasse2016; Rodgers et al.2008; Hauck and Völker2015). Over the 21st century and under a high-greenhouse-gas-emission scenario, Earth system model simulations project that the seasonal amplitude in surface [H+] will increase by 81 %, whereas the seasonal amplitude for aragonite saturation state (ΩA) is projected to decrease by 9 % globally on average (Kwiatkowski and Orr2018). Recent observation-based estimates as well as theoretical arguments support these projected increases in seasonality for [H+] and pCO2 (Landschützer et al.2018; Fassbender et al.2018). Thus, when extremes are defined with respect to a shifting baseline (i.e., mean state changes are neglected), the frequency and intensity of extreme [H+] events will likely increase due to increases in variability.

Unlike for marine heatwaves (Frölicher et al.2018; Collins et al.2020) and extreme sea level events (Oppenheimer et al.2020), little is known about the characteristics and changes of extreme ocean acidity events and, if so, only on seasonal timescales (Kwiatkowski and Orr2018). A global view of how extreme events in ocean chemistry will unfold in time and space and a mechanistic understanding of the relevant processes is missing. This knowledge gap is of particular concern as it is expected that extreme events in ocean acidity, defined with respect to both a fixed and a shifting baseline, are likely to become more frequent and intense under increasing atmospheric CO2. Given the potential for profound impacts on marine ecosystems, quantifying trends and patterns of extreme events in ocean acidity is a pressing issue.

In this study, we use daily mean output of a five-member ensemble simulation under low- and high-CO2-emission scenarios of a comprehensive Earth system model to investigate how the occurrence, intensity, duration, and volume of [H+] and Ω extreme events change under rising atmospheric CO2 levels. Extreme events defined with respect to both a fixed preindustrial and a shifting baseline are assessed, but the main focus is on extremes with respect to a shifting baseline and how these are affected by variability changes.

2 Methods

2.1 Model and experimental design

The simulations used in this study were made with the fully coupled carbon–climate Earth system model developed at the NOAA Geophysical Fluid Dynamics Laboratory (GFDL ESM2M) (Dunne et al.2012, 2013). The GFDL ESM2M model consists of ocean, atmosphere, sea-ice, and land modules and includes land and ocean biogeochemistry. The ocean component is the Modular Ocean Model version 4p1 (MOM4p1), with a nominal 1 horizontal resolution increasing to 1∕3 meridionally at the Equator, with a tripolar grid north of 65 N, and with 50 vertical depth levels. The MOM4p1 model has a free surface, with the surface level centered around about 5 m depth, and the spacing between consecutive levels is about 10 m down to a depth of about 230 m (Griffies2009) with increasing spacing below. The dynamical sea-ice model uses the same tripolar grid as MOM4p1 (Winton2000). The atmospheric model version 2 (AM2) has a horizontal resolution of 2×2.5 with 24 vertical levels (Anderson et al.2004). The land model version 3 (LM3) simulates the cycling of water, energy, and carbon dynamically and uses the same horizontal grid as AM2 (Shevliakova et al.2009).

The ocean biogeochemical and ecological component is version two of the Tracers of Ocean Phytoplankton with Allometric Zooplankton (TOPAZv2) module that parameterizes the cycling of carbon, nitrogen, phosphorus, silicon, iron, oxygen, alkalinity, lithogenic material, and surface sediment calcite (see Supplement in Dunne et al.2013). TOPAZv2 includes three explicit phytoplankton groups – small, large, and diazotrophs – and one implicit zooplankton group. The ocean carbonate chemistry is based on the OCMIP2 parameterizations (Najjar and Orr1998). The dissociation constants for carbonic acid and bicarbonate ions are from Dickson and Millero (1987), which are based on Mehrbach et al. (1973), and the carbon dioxide solubility is calculated according to Weiss (1974). Total alkalinity in TOPAZv2 includes contributions from phosphoric and silicic acids and their conjugate bases. TOPAZv2 also simulates diurnal variability in ocean physics as well as in phytoplankton growth. While diurnal variations in open-ocean pH are therefore simulated to some extent, we do not expect the model to fully capture the high diurnal variability in seawater chemistry, especially in coastal regions with large biological activity (Kwiatkowski et al.2016; Hofmann et al.2011), due to its relatively coarse resolution and simple biogeochemical model.

We ran a five-member ensemble simulation covering the historical 1861–2005 period, followed by a high (RCP8.5; RCP: Representative Concentration Pathway) and a low-greenhouse-gas-emission scenario (RCP2.6) over the 2006–2100 period with prescribed atmospheric CO2 concentrations. RCP8.5 is a high-emission scenario without effective climate policies, leading to continued and sustained growth in greenhouse gas emissions (Riahi et al.2011). In the GFDL ESM2M model, global atmospheric surface temperature in the RCP8.5 ensemble is projected to increase by 3.24 C (ensemble minimum of 3.17 C to ensemble maximum of 3.28 C) between the preindustrial period and 2081–2100. The RCP2.6 scenario represents a low-emission, high-mitigation future (van Vuuren et al.2011) with a simulated warming in the GFDL ESM2M model of 1.21 (1.18–1.26) C. The five ensemble members over the historical period were initialized from a multicentury preindustrial control simulation that was extended with historical land use over the 1700–1860 period (Sentman et al.2011). The five ensemble members were generated by adding different sea surface temperature (SST) disturbances of the order 10−5C to a surface grid cell in the Weddell Sea at 70.5 S, 51.5 W on 1 January 1861 (Wittenberg et al.2014; Palter et al.2018). Although the ocean biogeochemistry is not perturbed directly, [H+] and Ω differences between the ensemble members spread rapidly over the globe. On average, the ensemble members can be regarded as independent climate realizations after about 3 years of simulation for surface waters and after about 8 years at 200 m (Frölicher et al.2020). Neither the choice of the perturbation location nor the choice of the perturbed variable has a discernible effect on the results presented here (Wittenberg et al.2014). In addition, an accompanying 500-year preindustrial control simulation was performed.

2.2 Analysis

2.2.1 Extreme-event definition and characterization

We analyze daily mean data of [H+] and aragonite saturation state ΩA in the upper 200 m of the water column. [H+] is on the total scale and hence the sum of the concentrations of free protons and hydrogen sulfate ions. We define an event as a [H+] extreme event when the daily mean [H+] exceeds the 99th percentile, i.e., occurring once every 100 d. Similarly, we define a ΩA extreme event when the daily mean ΩA falls below the 1st percentile. The percentiles are calculated for each grid cell from daily mean data of the 500-year preindustrial control simulation. In contrast to absolute thresholds, relative thresholds, such as those used here, take into account regional differences in a variable's mean state, variance, and higher moments. Events that are defined based on relative thresholds have the same probability of occurrence across the globe in the period in which they are defined (e.g., preindustrial period; see also Frölicher et al.2018).

We assess changes in [H+] and ΩA extreme events when they are defined with respect to both a fixed preindustrial baseline and a shifting baseline. Under the fixed baseline approach, the secular trends as well as changes in variability and the higher moments of the distribution impose changes in extreme events. Under the shifting-baseline approach, which is the focus of this study, a value is considered extreme when it is much higher or lower than the baseline that undergoes changes due to secular trends in the variable. Thus, changes in the different extreme-event characteristics are only caused by changes in variability and the higher moments of the distributions. To define the extreme events with respect to the shifting baselines, we subtract the secular trends in [H+] and ΩA at each grid cell and in each individual ensemble member prior to the calculation of the different extreme-event characteristics based on the preindustrial percentiles (depicted for one grid cell in Fig. 1). The secular trend is calculated as the five-member ensemble mean, which has been additionally smoothed with a 365 d running mean to keep the seasonal signal in the data (further information in Appendix A). The removal of the secular trend ensures that the mean state in the processed data stays approximately constant while day-to-day to interannual variability can change over the simulation period (Fig. 1).

We calculate four extreme-event metrics: (a) the yearly extreme days (in days; number of days per year above the 99th percentile for [H+] and below the 1st percentile for ΩA), (b) the annual mean duration (in days; the average number of days above the 99th percentile for [H+] and below the 1st percentile for ΩA of single events within a year), (c) the annual mean maximal intensity (in nmol kg−1 or ΩA unit; maximum [H+] or ΩA anomalies with respect to the percentile threshold over the duration of a single extreme event and then averaged over all events within a year), and (d) the mean volume covered by individual extreme events in the upper 200 m (in km3; mean volume of 3D clusters of connected grid cells that are above the 99th percentile for [H+] or below the 1st percentile for ΩA, calculated using the measure.label function from the scikit-image library for Python for each day; these daily means are then averaged annually). The yearly extreme days, duration, and maximal intensity are calculated for individual grid cells at the surface and at 200 m. While the truncation of extremes between years alters the results for duration and maximal intensity, it allows for the calculation of annual extreme-event characteristics. We focus our analysis not only on the surface, but also on 200 m to study changes in extreme events within the seasonal thermocline. Most organisms susceptible to ocean acidification are found in the upper 200 m, such as reef-forming corals and calcifying phytoplankton.

2.2.2 Decomposition of [H+] variability into different variability components

We use three steps to decompose the total temporal variability in [H+] into interannual, seasonal, and subannual variability (Fig. 2). In a first step, we calculate the climatological seasonal cycle from the daily mean data by averaging each calendar day over all years in the time period of interest. Seasonal variability is then identified with the time-series variance of this 365 d long seasonal cycle. The secular trend in the daily mean data has been removed with the five-member ensemble mean before doing the analysis. In a second step, we subtract the seasonal cycle from the data and estimate the spectral density (Chatfield1996) of this residual time series using the periodogram function from the scipy.signal Python library. In a third step, we calculate the variance arising from variations on interannual and subannual timescales from the spectral density to obtain interannual and subannual variability (further information is given in Appendix B). Following this methodology, subannual variability comprises all variations in daily mean data with periodicities of less than a year that are not part of the seasonal cycle.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f01

Figure 1Simulated daily mean surface [H+] (a) and ΩA (c) at 40N and 30W in the North Atlantic for one ensemble member over the preindustrial period, the 1861–2005 historical period, and the 2006–2100 period under RCP8.5. The same data as in (a) and (c) but with subtracted ensemble-mean changes with respect to the average of the 500-year preindustrial control simulation is shown in panels (b) and (d). For [H+], the preindustrial 99th percentile threshold (horizontal blue line in panels a and b) is increasingly exceeded even when subtracting the ensemble-mean change, because [H+] variability increases. In contrast, a reduction in ΩA variability leads to a reduced undershooting of the preindustrial 1st percentile (d).

Download

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f02

Figure 2The three-step decomposition of [H+] variance into interannual, seasonal, and subannual variance, exemplified for a surface grid cell at 40 N and 30 W in the North Atlantic in the preindustrial control simulation. In a first step, the climatological seasonal cycle is determined (over the whole period, only 5 years are depicted here) and its variance is calculated. Note that the seasonal cycle in this grid cell has two minima and maxima. In a second step, the spectral density of the anomalies with respect to the seasonal cycle is calculated. In a third step, interannual and subannual variance is estimated from the spectral density.

Download

2.2.3 Taylor expansion of [H+] and ΩA variability changes

To understand the processes behind the simulated changes in the variabilities of [H+] and ΩA, we decompose these changes into contributions from changes in temperature (T), salinity (S), total alkalinity (AT), and total dissolved inorganic carbon (CT). Assuming linearity, the difference of [H+] from its mean at time step i can be decomposed into contributions from the drivers by employing a first-order Taylor expansion,

(1) H + ( i ) - H + H + C T C T , A T , T , S C T ( i ) - C T + H + A T C T , A T , T , S A T ( i ) - A T + H + T C T , A T , T , S T ( i ) - T + H + S C T , A T , T , S S ( i ) - S ,

and analogously for ΩA. The partial derivatives are evaluated at T, S, CT, and AT, i.e., the temporal mean values of the drivers in the period of interest. While it is important to take into account the climatological total phosphate and total silicate concentrations for calculating the partial derivatives (Orr and Epitalon2015), one introduces only small errors by neglecting variations in phosphate and silicate. The partial derivatives in Eq. (1) are evaluated using mocsy 2.0 (Orr and Epitalon2015).

Using the Taylor decomposition (Eq. 1), one can for example express the seasonal variation in [H+] as a function of the drivers' seasonal variations (Kwiatkowski and Orr2018). In this study, however, we analyze the time-series variance of [H+] and ΩA that also includes variability on other timescales (see Sect. 2.2.2) and the drivers of its changes. From the Taylor approximation (Eq. 1) and the definition of variance (e.g., Coles2001), it follows that the variance of [H+] can be written as a function of the partial derivatives with respect to the drivers (sensitivities), the standard deviations of the drivers, and their pairwise correlation coefficients:

(2) σ H + 2 = H + C T 2 σ C T 2 + H + A T 2 σ A T 2 + H + T 2 σ T 2 + H + S 2 σ S 2 + 2 H + C T H + A T cov C T , A T + 2 H + C T H + T cov C T , T + 2 H + C T H + S cov C T , S + 2 H + A T H + T cov A T , T + 2 H + A T H + S cov A T , S + 2 H + T H + S cov T , S ,

where the pairwise covariances are functions of the standard deviations and correlation coefficients according to cov(x,y)=σxσyρx,y, and the partial derivatives are again evaluated at the temporal mean values T, S, CT, and AT. This methodology has also been used to propagate uncertainties in carbonate system calculations (Dickson and Riley1978; Orr et al.2018) and to identify drivers of potential predictability in carbonate system variables (Frölicher et al.2020). Based on Eq. (2) and the analogous result for ΩA, a change in variance of [H+] and ΩA can be attributed to changes in the sensitivities that arise from changes in the drivers' mean states, to changes in the drivers' standard deviations, and to changes in the pairwise correlations between the drivers. We do so by calculating the Taylor series of Eq. (2) (further information in Appendix C). We then identify the [H+] variance change from mean changes in the drivers as the sum of all terms in the expansion that describe the contributions of sensitivity changes to the overall change in variance (ΔsσH+2). Likewise, we identify the contribution from standard deviation changes in the drivers (ΔσσH+2). We further group terms in the expansion that stem from simultaneous changes in the sensitivities and standard deviations (ΔsσσH+2) and the remaining terms that arise either from correlation changes alone or mixed contributions from correlation changes and changes in sensitivities and standard deviations (Δρ+σH+2). Since these four components contain all terms in the Taylor series, they exactly reproduce a change in variance represented by Eq. (2),

(3) Δ σ H + 2 = Δ s σ H + 2 + Δ σ σ H + 2 + Δ s σ σ H + 2 + Δ ρ + σ H + 2 .

We also assess the contributions to the four components from CT alone; from CT and AT; and from CT, AT, and T. The equivalent procedure is also used to decompose variance change in ΩA. Further information on the decomposition is given in Appendix C.

2.3 Model evaluation

The focus of our analysis is on changes in variability in [H+] and ΩA. As observation-based daily mean data of the inorganic carbon chemistry at the global scale are not available, we limit the evaluation of the Earth system model simulation to the representation of the seasonal cycles of [H+] and ΩA, and especially on its changes over the 1982–2015 period. We developed an observation-based dataset for surface monthly [H+] and ΩA using monthly surface salinity, temperature, pCO2, and AT fields. Salinity and temperature data are taken from the Hadley Centre EN.4.2.1 analysis product (Good et al.2013). AT is then calculated using the LIARv2 total alkalinity regression from salinity and temperature (Carter et al.2018). For pCO2, we use the neural-network-interpolated monthly data from Landschützer et al. (2016), which is based on SOCATv4 (Bakker et al.2016). Although not fully capturing pCO2 variability in regions with only few observations (Landschützer et al.2016), the pCO2 dataset appears to be generally well suited for analyzing pCO2 seasonality and changes therein (Landschützer et al.2018). An exception is the Southern Ocean, where data-based pCO2 products are uncertain due to sparse data in winter (Gray et al.2018). [H+] and ΩA are then calculated from salinity, temperature, AT, and pCO2 using the CO2SYS carbonate chemistry package (van Heuven et al.2011). Uncertainties in the derived seasonal cycles for [H+] and ΩA that arise from uncertainties in the observation-based input variables are not quantified in this study.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f03

Figure 3Seasonal amplitude of [H+] (calculated as yearly maximum minus the yearly minimum after subtracting a cubic spline from the data) over the period 1982–2015 averaged over five different latitude bands for the observation-based estimate (a) and the GFDL ESM2M model historical (1982–2005) and RCP8.5 (2006–2015) ensemble simulations (b), along with the same for data-based ΩA (c) and simulated ΩA (d). Linear trends in all panels are overlaid as thick lines. The linear trend of the simulated changes is calculated as the mean of the five individual ensemble trends.

Download

In most regions, the GFDL ESM2M model captures the observation-based mean seasonal cycle in [H+] and ΩA quite well, in particular for ΩA (the mean values of the seasonal amplitudes in Fig. 3). However, potential biases in the mean seasonal amplitudes do not directly have an effect on projected changes in extreme events, as we base the extreme-event definition on relative thresholds.

We then compare the simulated ensemble-mean trends in seasonal amplitude with the observation-based estimates (Fig. 3; Appendix D). Similar to the mean seasonal cycle results, the GFDL ESM2M model captures the observed trends in the seasonal [H+] and ΩA amplitudes for different latitudinal bands over the 1982–2015 period relatively well. The ensemble-mean trends in the simulated seasonal [H+] amplitudes are positive for all latitude bands (Fig. 3, Table 1), consistent with the observation-based estimates. While the estimates for the simulated trends are significantly larger than zero for all latitude bands, this is not the case for the observation-based trends in the equatorial region (10 S–10 N) and the northern low latitudes (10–40 N) (Table 1). The simulated [H+] seasonality trends are significantly smaller (with 90 % confidence level) than estimated from observations in the northern high (40–80 N; orange thick lines in Fig. 3a, b) and southern low latitudes (40–10 S; blue thick lines in Fig. 3a, b), where the trends from the model ensemble are 0.031±0.012 nmol kg−1 per decade and 0.035±0.003 nmol kg−1 per decade, compared to the observation-based trends of 0.106±0.040 nmol kg−1 per decade and 0.055±0.014 nmol kg−1 per decade, respectively. The simulated ensemble-mean trends for the remaining latitude bands are not significantly different from the observation-based trend estimates.

For the seasonal amplitude of ΩA, we find a significant negative trend in the observation-based data in the northern low latitudes and significant negative trends in the simulations in the northern and southern high latitudes (Table 1). The negative trends in seasonal amplitude in the simulations are significantly different from the observation-based trends in the northern high latitudes (-0.015±0.004 vs. 0.002±0.009 ΩA units per decade) and in the southern high latitudes (-0.012±0.002 vs. 0.000±0.005 ΩA units per decade).

Table 1Linear trends in seasonal amplitude of [H+] (in nmol kg−1 per decade) and ΩA (in 10−3 per decade) for five latitude bands over the period 1982–2015. Results are shown for the observation-based data (Obs.) and the five-member ensemble mean of the ESM2M model simulations (ESM2M) following the RCP8.5 scenario over 2006–2015. The range (±) denotes the 90 % confidence interval.

Download Print Version | Download XLSX

In summary, taking into account previous evaluations of the mean states of [H+] and ΩA and the underlying drivers in the GFDL-ESM2M model (Bopp et al.2013; Kwiatkowski and Orr2018), the model performs well against a number of key seasonal performance metrics. However, the model slightly underestimates past increases in seasonal amplitude of [H+], especially in the northern and southern high latitudes. In contrast to the observation-based data, the model also projects negative trends in the ΩA seasonal amplitude there. Nevertheless, the observation-based trends in the northern and especially southern high latitudes are uncertain because wintertime data are sparse there. Even though we lack the daily mean observation-based data to undertake a full assessment, it appears that the GFDL ESM2M model is adequate to assess changes in open-ocean variability of [H+] and ΩA and to assess changes in extreme events that arise thereof.

3 Results

We first briefly discuss the simulated changes in [H+] and ΩA extreme events when these events are defined with respect to a fixed preindustrial baseline period (Sect. 3.1). In Sect. 3.2 and 3.3, these results are contrasted with changes in extremes that are defined with respect to a shifting baseline, i.e., where the secular trends do not alter extreme events. In Sect. 3.4, variability changes are decomposed into subannual, seasonal, and interannual variability contributions. The processes leading to variability changes are analyzed in Sect. 3.5.

3.1 Global changes in extremes defined relative to a fixed preindustrial baseline

When using the fixed preindustrial 99th and 1st percentiles to define extreme events in [H+] and ΩA, respectively, large increases in the number of days with [H+] and ΩA extremes are projected over the 1861–2100 period in both low- and high-CO2-emission scenarios (Figs. 4 and A1). Over the historical period, the GFDL ESM2M model projects an increase in yearly extreme days for surface [H+] from 3.65 d per year during the preindustrial period to 299 d per year in 1986–2005. By year 2030 and under both CO2 emission scenarios, the surface ocean is projected to experience a “near-permanent acidity extreme state”; i.e., [H+] is above the preindustrial 99th percentile more than 360 d pear year. Likewise, the average duration of events saturates near 365 d, and the intensity of events increases strongly, mainly reflecting the large increase in mean [H+] (Fig. A1). A similar but slightly delayed evolution in the number, maximal intensity, and duration of [H+] extremes is simulated at 200 m (Fig. A1).

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f04

Figure 4Simulated globally averaged yearly extreme days defined with respect to a fixed baseline for [H+] using the preindustrial 99th percentile (a) and for ΩA using the preindustrial 1st percentile (b). Shown are changes at the surface over the 1881–2100 period following historical (black lines) and future scenarios, RCP8.5 (red) and RCP2.6 (blue). The thick lines display the five-member ensemble means, and the shaded areas represent the maximum and minimum ranges of the individual ensemble members.

Download

Large increases in yearly extreme days are also projected for ΩA when using a fixed preindustrial 1st percentile as a baseline (Fig. 4b). Similar to [H+], the entire surface ocean is projected to approach a permanent ΩA extreme state during the 21st century under the RCP8.5 scenario. A near-permanent extreme state is projected by year 2062. In contrast to [H+], a permanent ΩA extreme state of the global ocean is avoided under the RCP2.6 scenario.

3.2 Global changes in extremes defined relative to a shifting baseline

Next, we investigate changes in [H+] and ΩA extremes when the extreme events are defined with respect to a shifting (time-moving) baseline; i.e., changes in extremes arise only from changes in variability and higher moments of the distributions. The GFDL ESM2M model projects large increases in the number, intensity, duration, and volume of [H+] extreme events over the 1861–2100 period (Fig. 5). Over the historical period (from the preindustrial period to 1986–2005), the model projects that the number of surface [H+] extreme days increases from 3.65 d per year to 10.0 d per year (Fig. 5a, ensemble ranges are given in Table 2). The maximal intensity is projected to increase from 0.08 to 0.12 nmol kg−1 (Fig. 5c, Table 2) and the duration from 11 to 15 d (Fig. 5e). Compared to preindustrial conditions, these changes correspond to a 173 % increase in the number of days per year, a 44 % increase in the maximal intensity, and a 45 % increase in the duration of [H+] extreme events. The volume of individual events is projected to increase by 20 % over the historical period, from a typical volume of 2.7×103 km3, which is about 0.004 % of the total ocean volume in the upper 200 m (Fig. 5g), to 3.2×103 km3.

Over the 21st century, extreme events in ocean acidity, defined with respect to a shifting baseline, are projected to further increase in frequency, intensity, duration, and volume (Fig. 5). By 2081–2100 under the RCP8.5 scenario, the number of [H+] extreme days per year at the surface is projected to increase to 50 d (corresponding to a 1273 % increase relative to the preindustrial period). The maximal intensity is projected to increase to 0.38 nmol kg−1 (371 % increase), the duration to 32 d (199 % increase), and the volume to 13.9×103 km3 (414 % increase).

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f05

Figure 5Simulated changes in globally averaged [H+] extreme-event characteristics over the 1881–2100 period following historical (black lines) and future RCP8.5 (red) and RCP2.6 (blue) scenarios. The extreme events are defined with respect to a shifting baseline. Yearly extreme days, maximal intensity, and duration are shown for the surface (a, c, e) and for 200 m (b, d, f). Volume is shown in (g). The thick lines display the five-member ensemble means, and the shaded areas represent the maximum and minimum ranges of the individual ensemble members.

Download

At 200 m, the [H+] extreme events in preindustrial conditions are in general more intense (0.17 nmol kg−1; Fig. 5d) and longer lasting (38 d; Fig. 5f) than at the surface. The stronger extreme events are caused by the overall larger variability at 200 m than at the surface in the preindustrial period. The longer duration is connected to the more pronounced contribution from interannual variability (see Sect. 3.4). However, projected relative changes over the historical period and the 21st century are smaller at 200 m than at the surface and with larger year-to-year variations across the ensembles. Under recent past conditions (1986–2005), the number of extreme days per year at 200 m is 4.3 d per year (corresponding to an 18 % increase since the preindustrial period), the maximal intensity 0.20 nmol kg−1 (18 % increase), and the duration 46 d (21 % increase). By the end of the 21st century under the RCP8.5 scenario, the number of [H+] extreme days per year is projected to increase to 32 d per year, the maximal intensity to 0.34 nmol kg−1, and the duration to 99 d. Notably, extreme events in [H+] are projected to become slightly less intense at 200 m than at the surface (0.34 vs. 0.38 nmol kg−1) by the end of the century under RCP8.5, even though they were more intense in preindustrial times at depth. In contrast, surface [H+] extreme events remain shorter in duration at the end of the century than at 200 m.

Table 2Simulated global ensemble-mean [H+] extreme-event characteristics, when extremes are defined with respect to a shifting baseline. Values in brackets denote ensemble minima and maxima.

Download Print Version | Download XLSX

Under the RCP2.6 scenario, the magnitude of changes in the different [H+] extreme-event characteristics by the end of the century is substantially smaller than in the RCP8.5 scenario. This difference is especially pronounced at the surface (blue lines in Fig. 5). There, the number of extreme days per year, maximal intensity, and duration under the RCP2.6 are projected to be 46 % (44–47), 43 % (43–44), and 75 % (73–77) of that under the RCP8.5 scenario. At depth, the differences between the RCP2.6 and RCP8.5 scenario are less pronounced and only emerge in the second half of the 21st century. As opposed to the surface, the number of [H+] extreme days per year and the maximal intensity at 200 m as well as the volume of events are projected to increase significantly even after the atmospheric CO2 concentration stabilizes in RCP2.6 around year 2050. This delayed response in the subsurface is due to the relatively slow surface-to-subsurface transport of carbon. However, this is not the case for the duration, which slightly decreases in the second half of the 21st century at depth (Fig. 5f). This decrease in duration mainly occurs in the subtropics, where events generally last longer (Fig. A3b). It is connected to an increase in the contribution from high-frequency variability to total variability in those regions over that period.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f06

Figure 6Simulated changes in the yearly number of ΩA extreme days. The extreme events are defined with respect to a shifting baseline. Panels (a) and (b) show the globally averaged simulated yearly extreme days in ΩA from 1881 to 2100 following historical (black lines) and future RCP2.6 (blue) and RCP8.5 (red) scenarios at the surface (a) and 200 m (b). The thick lines display the five-member ensemble means, and the shaded areas represent the maximum and minimum range of the individual ensemble members. Panels (c) and (d) show the simulated regional changes in yearly extreme days in ΩA from the preindustrial period to 2081–2100 under the RCP8.5 scenario at the surface (c) and at 200 m (d). Shown are changes averaged over all five ensemble members. The black lines highlight the pattern structure and gray colors represent regions where no ensemble member simulates extremes during 2081–2100.

In contrast to [H+] extreme events, the number of yearly extreme days in ΩA is projected to decrease over the historical period and during the 21st century under both the RCP8.5 and RCP2.6 scenarios (Fig. 6a–b, Appendix Table A1) when the extreme events are defined with respect to a shifting baseline. The number of surface ΩA extreme days per year by the end of this century is projected to be 63 % smaller under RCP8.5 and 39 % smaller under RCP2.6 compared to preindustrial conditions (ensemble ranges are given in Table A1). Projected changes at depth are less pronounced than at the surface, again with larger decreases under RCP8.5 than under RCP2.6. Despite this decline in extreme events when defined with respect to a shifting baseline, the long-term decline in the mean state of ΩA still leads to more frequent occurrence of extreme low ΩA events when defined with respect to a fixed baseline (see Sect. 3.1).

3.3 Regional changes in extremes defined relative to a shifting baseline

Surface [H+] extremes that are defined with respect to a shifting baseline are projected to become more frequent in 87 % of the surface ocean area by the end of the 21st century under the RCP8.5 scenario. However, the projected changes in these ocean acidity extremes are not uniform over the globe (Fig. 7; Appendix Fig. A3). The largest increases in the number of [H+] extreme days per year are projected in the Arctic Ocean (up to +120 d per year), in the subtropical gyres (up to +60 d per year), and in parts of the Southern Ocean and near Antarctica. There are also some regions including the eastern equatorial Pacific and parts of the Southern Ocean where the number of yearly extreme days in surface [H+] is projected to decrease. These are in general also the regions where the seasonality in [H+] is projected to decrease (see Sect. 3.4 below). The largest absolute changes in intensity of surface [H+] extremes (Fig. 7c) are projected for the subtropics, especially in the Northern Hemisphere. For example, events become up to 1 nmol kg−1 more intense in the subtropical North Pacific and Atlantic, corresponding roughly to a 10-fold increase in intensity with respect to the preindustrial period. The largest relative increases in intensity are projected for the Arctic Ocean, the North Atlantic, and around Antarctica, where more-than-10-fold increases with respect to the preindustrial period are projected. Regions with large increases in the number of yearly extreme days tend to also show large increases in the duration of extreme events (Fig. 7e). The Arctic Ocean is an exception. Although the number of yearly extreme days increases strongly, the increase in duration is not as pronounced. This discrepancy is because extremes are already long lasting but rare at preindustrial times (Fig. A3). So even though extreme events are projected to occur each year by the end of the century under RCP8.5, the increase in duration is relatively small.

At 200 m, the projected pattern of changes in yearly extreme days generally resembles that at the surface (Fig. 7b). The largest increases in yearly extreme days are projected for parts of the subtropics, the Southern Ocean, and the Arctic Ocean. In contrast to the surface, [H+] extremes at 200 m are projected to become less frequent in the equatorial Atlantic, the northern Indian Ocean, the North Pacific, and in large parts of the Southern Ocean. The regions indicating a decline in [H+] extremes at depth include also some of the eastern boundary current systems, such as the Humboldt, California, and Benguela Current systems. In most of these regions, extreme events are projected to disappear in the RCP8.5 scenario by the end of this century (gray regions in Fig. 7b). The largest increases in subsurface event intensity are projected in the subtropics (Fig. 7d), whereas the duration of [H+] extremes is projected to increase strongly in many regions of the mid-to-high latitudes of both hemispheres (Fig. 7f). The projected increases in duration at 200 m are much larger than at the surface.

The increase in the number of extreme days per year, the maximal intensity, and the duration is smaller under RCP2.6 compared to RCP8.5 for most of the ocean (Fig. A2). The largest increases in occurrence of extremes under RCP2.6 are simulated for the Arctic Ocean, similar to under RCP8.5, and for parts of the Southern Ocean. The regions in the Southern Ocean where the occurrence of extreme events in [H+] is projected to decrease largely overlap with those for RCP8.5, at the surface and at depth. On the other hand, unlike under RCP8.5, a decrease in extreme-event occurrence is only projected for a small fraction of the tropical oceans under RCP2.6.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f07

Figure 7Simulated regional changes in [H+] extreme-event characteristics from the preindustrial period to the 2081–2100 period under the RCP8.5 scenario at the surface and at depth for the yearly extreme days (a, b), the maximal intensity of events (c, d), and the duration of events (e, f). The extreme events are defined with respect to a shifting baseline. Shown are changes averaged over all five ensemble members. Gray colors represent areas where no extremes occur during 2081–2100, and the black lines highlight pattern structures.

While the decline in mean ΩA generally leads to lower values in ΩA and thus extreme events are becoming more frequent when defined with respect to a fixed preindustrial baseline (Sect. 3.1), extreme events in ΩA are projected to become less frequent throughout most of the ocean when defined with respect to a shifting baseline (89 % of surface area under RCP8.5 at the end of the 21st century; Fig. 6c). In many regions, extreme events in ΩA are projected to disappear by 2081–2100 under the RCP8.5 scenario (gray regions in Fig. 6c) when defined with respect to a shifting baseline. However, the number of yearly extreme days in ΩA is projected to increase by 10 or more in the subtropical gyres, especially in the western parts of the subtropical gyres. At 200 m, no extreme events are projected for most of the ocean during 2081–2100 under RCP8.5 (Fig. 6d).

3.4 Decomposition of temporal variability in [H+]

The changes in [H+] extreme events defined with respect to a shifting baseline mainly result from changes in [H+] variability. These variability changes may arise from changes in interannual variability, seasonal variability, and subannual variability. Thus, we decomposed the total variability into these three components (see Sect. 2.2.2). For the preindustrial period, the model simulates generally larger [H+] variance at depth than at the surface (0.42 vs. 0.15 nmol2 kg−2, not shown). Seasonality has the largest contribution at the surface (81 % of total variance). At 200 m, interannual variability has the largest contribution (63 %), and also subannual variability is more important compared to the surface (15 % vs. 8 %).

Over the 1861–2100 period under the historical-RCP8.5 forcing, changes in seasonality clearly dominate the overall change in variability at the surface with 87 % contribution to the overall variance change in the global mean (Fig. 8b, d). Changes in interannual variability (3 % contribution to overall variance change; Fig. 8a, d) and subannual variability (10 %; Fig. 8c, d) play a minor role. The largest increases in variability for all three variability types are projected for the northern high latitudes. Around Antarctica and the southern end of South America, large increases in seasonal variability are projected (Fig. 8b). In the tropical Pacific and parts of the Southern Ocean, decreases in interannual and seasonal variability are projected (Fig. 8a, b).

In contrast to the surface, changes in interannual and to a lesser extent subannual variability at 200 m are also important for explaining the overall changes in [H+] variability (Fig. 8e, g, h). Changes in interannual variability contribute most to overall variance change at the global scale (with 42 % contribution). Seasonal variability changes are almost equally important (37 %), and changes in subannual variability also contribute substantially to changes in total variability (20 %). The patterns of variability changes are very similar across the three temporal components of variability. The largest increases in [H+] variability are simulated north and south of the Equator. These regions tend to be already more variable during the preindustrial period (see Fig. A3a). However, the model also projects an increase in variability for regions that are less variable during the preindustrial period, such as northern high latitudes. All three temporal components of variability are projected to decrease in the tropics and parts of the Southern Ocean. The variability decrease in those regions is most pronounced for interannual variability (Fig. 8e).

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f08

Figure 8Contribution to projected changes in [H+] variance from interannual variability (a, e), seasonal variability (b, f), and subannual variability (c, g) between the preindustrial period and the 2081–2100 period under the RCP8.5 forcing at the surface (a–d) and at 200 m (e–h). Shown are the ensemble-mean changes. The black lines highlight the pattern structure. Zonal mean contributions are shown for the surface (d) and for 200 m (h). The sum of the three components (black lines) accurately reproduces the simulated variance change (gray dashed lines).

3.5 Drivers of [H+] and ΩA variability changes

In this section, we investigate the changes in the drivers that cause the variability changes in [H+] and ΩA. Drivers are carbon (CT), alkalinity (AT), temperature, and salinity. To do so, we attribute changes in [H+] and ΩA variability to four factors (see Sect. 2.2.3 for further details): (i) changes in the mean states of the drivers that control the sensitivities (ΔsσH+2), (ii) changes in the variabilities of the drivers (ΔσσH+2), (iii) simultaneous changes in the mean states and variabilities of the drivers (ΔsσσH+2; this contribution arises because both mean states and variabilities change and cannot be attributed to either (i) or (ii) alone), and (iv) changes in the correlations between the drivers, also including mixed contributions from correlation changes together with mean state and variability changes (Δρ+σH+2). In other words, (iv) describes the change in variability that arises because the correlations between the drivers also change, and not only their mean states and variabilities.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f09

Figure 9Decomposition of surface [H+] variability changes into different drivers (CT, AT, temperature, and salinity). Shown are changes from the preindustrial period to 2081–2100 following the RCP8.5 scenario. The simulated change in [H+] variance (ΔσH+2) is decomposed into the contribution from changes in the sensitivities that arise from changes in the drivers' mean values (ΔsσH+2), the contribution from changes in the drivers' standard deviations (ΔσσH+2), the contribution from simultaneous changes in the sensitivities and the drivers' standard deviations (ΔsσσH+2), and the contribution from correlation changes alone together with simultaneous changes in correlations and sensitivities and standard deviations (Δρ+σH+2) (a). The small mismatch between the sum of the components (black line) and simulated variance change (gray dashed line) arises because the decomposition is based on Eq. (2), which is an approximation to simulated [H+] variance. The contributions to these components from changes in CT alone (light blue lines); from changes in CT and AT (green lines); and from CT, AT, and temperature (gold lines) are shown in panels (b)(e). The dashed black lines in panels (b)(e) show the total components that contain contributions from all four drivers.

Download

The drivers' mean changes between the preindustrial period and 2081–2100 under RCP8.5 cause a strong increase in surface [H+] variability, which is most pronounced in the high latitudes (ΔsσH+2; red line in Fig. 9a, black dashed line in Fig. 9b). On a global average, these variance changes due to the mean changes in the drivers (ΔsσH+2=1.3 nmol2 kg−2) are much larger than the total simulated variance change in [H+] (ΔσH+2=0.5 nmol2 kg−2, dashed gray or solid black line in Fig. 9a). In general, an increase in mean CT, temperature, and salinity would lead to an increase in ΔsσH+2, whereas an increase in mean AT would lead to a decrease. The GFDL ESM2M model projects an increase in mean CT over the entire surface ocean (Fig. A5a) due to the uptake of anthropogenic CO2 from the atmosphere and therefore an increase in ΔsσH+2 (light blue line in Fig. 9b). In the high latitudes, a relatively small increase in mean CT leads to a large increase in ΔsσH+2, because [H+] is more sensitive to changes in CT due to the low buffer capacity there. Decreases in mean AT further contribute to the increase in ΔsσH+2 in the high latitudes (green line in Fig. 9b). In the low-to-mid latitudes and in particular in the Atlantic Ocean, mean surface AT is projected to increase (Fig. A5) and therefore dampens the overall increase in ΔsσH+2 (green line in Fig. 9b). The changes in AT are largely due to changes in freshwater cycling that also manifest in salinity changes (Fig. A5, Carter et al.2016). Increases in temperature additionally increase ΔsσH+2, mainly in the northern mid-to-high latitudes (gold line in Fig. 9b), but the overall impact of mean changes in temperature, and especially salinity, is small.

Why is the increase in ΔσH+2 (gray dashed or black solid line in Fig. 9a) smaller than the increase from the mean changes (i.e., ΔsσH+2; red line in Fig. 9a)? In the high latitudes, the projected change in the variability of the drivers (Fig. A6) contributes negatively to the [H+] variability change and counteracts to some degree the increase in ΔsσH+2. These variability changes alone have a small imprint on ΔσσH+2 (blue line in Fig. 9a; black dashed line in Fig. 9c), but the variability changes dampen the increases from the mean changes (ΔsσσH+2, magenta line in Fig. 9a, black dashed line in Fig. 9d). The latter contribution is large in the high latitudes, where mean changes alone would lead to a strong increase. In the high latitudes, decreases in CT variability (Fig. A6a) together with increases in mean CT (Fig. A5a) can explain much of the negative contribution from ΔsσσH+2 (light blue line in Fig. 9d). In the northern high latitudes, mean and variability changes in AT are also important for ΔsσσH+2 (green line in Fig. 9d). The additional contribution from changes in the correlations between the drivers (Δρ+σH+2; cyan line in Fig. 9a) also tends to contribute negatively to [H+] variability changes, especially in the North Atlantic, and changes in correlations with temperature play an important role (gold line in Fig. 9e). In summary, the increase in [H+] variability at the surface is mainly caused by increases in mean CT attenuated by decreases in CT variability in the high latitudes. Mean changes in AT reinforce the increase in [H+] variability in the northern high latitudes but dampen the increase in the low latitudes.

At 200 m, the projected increase in ΔσH+2 (gray dashed or black solid line in Fig. 10a) is also a result of the large increase due to the mean changes in the drivers (ΔsσH+2; red line in Fig. 10a; dashed black line in Fig. 10b) and the decrease due to the interplay between mean changes and decreases in the variability (ΔsσσH+2; magenta line in Fig. 10a, black dashed line in Fig. 10d). Similar to the surface, the changes in mean and variability of CT are the most important drivers of changes (light blue lines in Fig. 10b, d). Increases in mean AT partially compensate for the increase in [H+] variability due to the increase in mean CT (green lines in Fig. 10b, d). Changes in [H+] variability due to changes in temperature and salinity are small. In contrast to the surface, the individual compensating contributions to [H+] variability change from mean changes and simultaneous mean and variability changes in the drivers, in particular those in CT, are much larger at 200 m. The global average variance change due to the mean changes in the drivers (ΔsσH+2=3.7 nmol2 kg−2) is much larger than the overall simulated variance change (ΔσH+2=0.1 nmol2 kg−2). The contribution from changes in the correlations between the drivers is overall small (cyan line in Fig. 10a) and stems mainly from changes in the correlation between CT and AT (Fig. 10e). Taken together, the increase in [H+] variability at 200 m mainly arises from the balance between increases in mean CT and decreases in CT variability. Increases in mean AT dampen these changes.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f10

Figure 10Decomposition of [H+] variability changes at 200 m into different drivers (CT, AT, temperature, and salinity). Shown are changes from the preindustrial period to 2081–2100 following the RCP8.5 scenario. The simulated change in [H+] variance (ΔσH+2) is decomposed into the contribution from changes in the sensitivities that arise from changes in the drivers' mean values (ΔsσH+2), the contribution from changes in the drivers' standard deviations (ΔσσH+2), the contribution from simultaneous changes in the sensitivities and the drivers' standard deviations (ΔsσσH+2), and the contribution from correlation changes alone together with simultaneous changes in correlations and sensitivities and standard deviations (Δρ+σH+2) (a). The contributions to these components from changes in CT alone (light blue lines); from changes in CT and AT (green lines); and from CT, AT, and temperature (gold lines) are shown in panels (b)(e). The dashed black lines in panels (b)(e) show the total components that contain contributions from all four drivers.

Download

Unlike for [H+], both mean changes (ΔsσΩ2; red lines in Fig. 11) and variability changes in the drivers (ΔσσΩ2: blue lines in Fig. 11) lead to a decrease in ΩA variability (ΔσΩ2; black dashed lines in Fig. 11). At 200 m, variability changes are even the dominant driver for reductions in ΩA variability. Simultaneous changes in means and variabilities (ΔsσσΩ2; purple lines in Fig. 11) contribute positively and dampen the reduction in ΩA variability from mean and variability changes alone. Mean and variability changes in CT are the main drivers for changes in ΩA variability as indicated by the tight relation between the dashed and solid red, blue, and purple lines in Fig. 11, in particular at 200 m. An exception is the northern high latitudes, where AT changes also play a substantial role at the surface (not shown). Correlation changes in the drivers (Δρ+σΩ2; cyan lines in Fig. 11) are of secondary importance and have the largest imprint in the northern mid-to-high latitudes at the surface.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f11

Figure 11Decomposition of ΩA variability changes into different drivers. The simulated zonal mean contribution to variance changes in ΩA (black dashed lines, ΔσΩ2) from the preindustrial period to 2081–2100 (RCP8.5) at the surface (a) and at 200 m (b). Shown is the contribution from sensitivity changes (due to mean changes in the drivers) (red lines, ΔsσΩ2), standard deviation changes in the drivers (blue lines, ΔσσΩ2), simultaneous changes in sensitivities and standard deviations (purple lines, ΔsσσΩ2), and all contributions that involve changes in the drivers' correlations (cyan lines, Δρ+σΩ2). Furthermore, contributions from mean changes, standard deviation changes, and simultaneous mean and standard deviation changes in CT alone are shown (dashed red, blue, and purple lines, respectively).

Download

4 Discussion and conclusions

We provide a first quantification of the historical and future changes in extreme events in ocean acidity by analyzing daily mean three-dimensional output from a five-member ensemble simulation of a comprehensive Earth system model. In our analysis, we focus on changes in high-[H+] and low-ΩA extreme events that are defined with respect to a shifting baseline, where changes in extremes arise from changes in daily to interannual variability. Secular trends in the mean state were removed from the model output before analyzing extremes under this approach. We show that such extreme events in [H+] are projected to become more frequent, longer lasting, more intense, and spatially more extensive under increasing atmospheric CO2 concentration, both at the surface and also within the thermocline. Under RCP2.6, the increase in these extreme-event characteristics is substantially smaller than under RCP8.5. The increase in [H+] variability is a consequence of increased sensitivity of [H+] to variations in its drivers. It is mainly driven by the projected increase in mean CT and additionally altered by changes in CT variability and AT mean and variability as well as changes in the correlations between the drivers. In contrast to [H+], variability of ΩA is projected to decline in the future. Therefore, extreme events in ΩA are projected to become less frequent in the future when defined with respect to a shifting baseline. The reason for the decline in variability is that ΩA, unlike [H+], becomes less sensitive to variations in the drivers with the mean increase in CT. Furthermore, the projected reductions in the drivers' variabilities, mainly in CT, further reduce ΩA variability.

The analysis of extreme events defined with respect to fixed preindustrial percentiles reveals that the secular trends in [H+] and ΩA are so large that they lead to year-round or almost-year-round extreme events in the upper 200 m over the entire globe by the end of the 21st century, even under the low-emission scenario RCP2.6. Extreme events are no longer temporally and spatially bounded events that arise due to the chaotic nature of the climate system but describe a permanent new state. Under the fixed baseline approach, the relative contribution of changes in variability or higher moments of the distribution to the changes in the number of extremes is small. For example, the number of yearly extreme days for surface [H+] over the 1986–2005 period under the shifting-baseline approach is only 3.8 % of that when defining the extreme events with respect to a fixed preindustrial baseline. This fraction differs regionally, reaching more than 10 % in the North Pacific, the North Atlantic, and the Arctic Ocean. However, we recall here that the changes in the number of [H+] extremes when defined with respect to a shifting baseline are large. These changes in variability may need to be taken into account when assessing the impacts of ocean acidity changes on marine organisms, especially when organisms are likely to adapt to the long-term mean changes but not to changes in variability.

We use the 99th percentile of the distribution from a preindustrial simulation for the definition of the extreme [H+] events (i.e., a one-in-a-hundred-days event at preindustrial levels), but the results may depend on the choice of this threshold. We tested the sensitivity of our results under the shifting-baseline approach by using also the 99.99th percentile threshold (i.e., a one-day-in-27.4-years event at preindustrial levels). The relative increase in the numbers of extreme [H+] days per year is larger for these rare extremes (Fig. 12). For example, nearly every second day with [H+] exceeding the 99th percentile (red solid lines in Fig. 12) is also a day with [H+] exceeding the 99.99th percentile (red dotted lines in Fig. 12) by the end of the 21st century under RCP8.5, both at the surface and at depth. As a result of this large relative increase in rare extremes, the model projects as many days with [H+] exceeding the 99.99th percentile by the end of the century under RCP8.5 (red dotted lines in Fig. 12) as it projects days exceeding the 99th percentile under RCP2.6 (blue solid lines in Fig. 12).

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f12

Figure 12Globally averaged number of yearly extreme days for [H+] over the historical (black lines), RCP2.6 (blue), and RCP8.5 (red) simulations for the preindustrial 99th (solid lines) and 99.99th percentile (dotted lines) at the surface (a) and 200 m (b). The extreme events are defined with respect to shifting baselines.

Download

The projected increase in [H+] variability and decrease in ΩA variability also alters the occurrence of extreme events based on absolute thresholds. An often used threshold is ΩA=1, below which seawater is corrosive with respect to the calcium carbonate mineral aragonite (Morse and Mackenzie1990). We assess the influence of the general decline in ΩA variability at the time where a grid cell falls below ΩA=1 for the first time. To do so, we compare these times within the historical and RCP8.5 ensemble to those for the hypothetical case where ΩA variability stays at the preindustrial level but mean ΩA undergoes the ensemble-mean evolution. We find that the decline in ΩA variability, which is observed in the historical and RCP8.5 ensemble, leads to an average delay of the first occurrence of undersaturation by about 11 years at the surface and about 16 years at 200 m. At the surface, these delays of undersaturation occur throughout the high latitudes (Fig. 13a). At depth, the delays are most pronounced in the tropics (Fig. 13b), but delays also occur in the high latitudes. Assuming unchanged seasonality, McNeil and Matear (2008) found that seasonal aragonite undersaturation of surface waters in the Southern Ocean may occur 30 years earlier than annual mean aragonite undersaturation. However, our simulation shows that the reduction in ΩA variability delays the onset of undersaturation by about 10 to 15 years in the Southern Ocean relative to a hypothetical simulation where variability does not change. Therefore, changes in variability need to be taken into account when projecting the onset of seasonal undersaturation, especially in the high latitudes and in the thermocline of the tropics.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f13

Figure 13The difference in years between the first occurrence of aragonite undersaturation in the historical and RCP8.5 ensemble and a hypothetical simulation where variability does not change over the 1861–2100 period, but only the mean changes. Positive values (yellow and red) indicate a delayed onset of undersaturation resulting from declines in ΩA variability.

Previous studies have shown that the seasonal cycle of surface ocean pCO2 will be strongly amplified under increasing atmospheric CO2 (Gallego et al.2018; Landschützer et al.2018; McNeil and Sasse2016) and that a similar amplification is expected for surface [H+] (Kwiatkowski and Orr2018). Here we show that the changes in the seasonal cycle of [H+] translate into large increases in short-term extreme acidity events at the surface as well as at 200 m, when these events are defined with respect to a shifting baseline. In addition to earlier studies, we also show that changes in subannual variability, which are only partially resolved by monthly mean data, contribute to changes in extreme events in [H+] under increasing atmospheric CO2. Furthermore, we show that the average duration of extreme events at the surface and in recent past conditions (1986–2005) is about 15 d. To resolve such events that last for days to weeks, it is necessary to use daily mean output. Currently, ocean carbonate system variables from models that participate in the sixth phase of the Coupled Model Intercomparison Project are routinely stored with a monthly frequency on the Earth system grid (Jones et al.2016). We therefore recommend storing and using high-frequency output to study extreme events in the ocean carbonate systems.

Even though we consider our results as robust, a number of potential caveats remain. First, the horizontal resolution of the ocean model in the GFDL ESM2M model is rather coarse and cannot represent critical scales of small-scale circulation structures (e.g., Turi et al.2018). In addition, the biogeochemical processes included in the GFDL ESM2M model are designed for the open ocean but do not capture the highly variable coastal processes (Hofmann et al.2011). High-resolution ocean models with improved process representations are therefore needed to explore variability in ocean carbonate chemistry, especially in coastal regions and smaller ocean basins, such as the Arctic (Terhaar et al.2019a, b). Observation-based carbonate system data with daily mean resolution would also be necessary to thoroughly evaluate the models' capability to represent day-to-day variations in carbonate chemistry. Secondly, our results, in particular at the local scale, might depend on the model formulation. As the mean increases in CT mainly drive the increases in [H+] variability (see Fig. 9b), we expect that models with larger oceanic uptake of anthropogenic carbon show larger increases in [H+] variability than models with lower anthropogenic carbon uptake. The GFDL ESM2M model matches observation-based estimates of historical global anthropogenic CO2 uptake relatively well but still has difficulties in representing the regional patterns in storage (Frölicher et al.2015). Therefore, the exact regional patterns of CT changes may differ from model to model. Further studies focusing on the physical processes that lead to the regional CT changes may help to better constrain the regional patterns in variability changes. In addition, it is currently rather uncertain how [H+] and ΩA variability changes as a result of changes in the drivers' variabilities. We have demonstrated that this factor is particularly important for ΩA and for [H+] at depth. It is well known that current Earth system models have imperfect or uncertain representations of ocean variability over a range of timescales (Keller et al.2014; Resplandy et al.2015; Frölicher et al.2016). A possible way forward would be to assess variability changes and changes in ocean acidity extreme events within a multimodel ensemble, which would likely provide upper and lower bounds. Finally, it is assumed that physical and biogeochemical changes in the ocean will also increase diurnal variability. In particular in coastal areas, such diurnal variations can have amplitudes that are much larger than the projected changes over the 21st century (Hofmann et al.2011). However, the GFDL ESM2M model does not fully resolve the diurnal variability. Future studies with Earth system models that resolve diurnal processes are needed to quantify changes in diurnal variability and the impacts of these changes on extreme acidity events.

Our results may also have important consequences for our understanding of the impacts of ocean acidification on marine organisms and ecosystems. The projected increase in the frequency and the duration of ocean acidity extremes implies that marine organisms will have less time to recover from high-[H+] events in the future. Organisms that cannot adapt to the large long-term changes in mean [H+] will likely be the most impacted. However, even if organisms may be able to adapt to the long-term increase in [H+], the large projected increase in [H+] extreme events due to changes in variability may push organisms and ecosystems to the limits of their resilience, especially those organisms that are commonly accustomed to a more steady environment. The risks for substantial ecosystem impacts are aggravated by the fact that the frequency and intensity of marine heatwaves are also projected to substantially increase (Frölicher et al.2018), which also negatively impact marine ecosystems (Wernberg et al.2016; Smale et al.2019). The interactions of intensified multiple stressors have the potential to influence marine ecosystems and the ocean’s biogeochemical cycles in an unprecedented manner (Gruber2011). However, further research is needed to understand the combined impacts of short-term ocean acidity extremes and marine heatwaves on marine ecosystems.

In conclusion, our analysis shows that [H+] and ΩA in the upper 200 m are projected to be almost permanently under extreme conditions by the end of the 21st century when extremes are defined relative to preindustrial baselines. Even when accounting for the changes in the long-term mean, short-term extreme events in [H+] are projected to become more frequent, to last longer, to be more intense, and to cover larger volumes of seawater due to increases in [H+] variability, potentially adding to the stress on organisms and ecosystems from the long-term increase in ocean acidity.

Appendix A: Identifying and removing the secular trend in the model data

In Sect. 3.2 and 3.3, we analyze the changes in extreme events in [H+] and ΩA that arise from day-to-day to interannual variability changes in these variables. We therefore need to remove the secular trends from the data prior to analysis. We estimate the secular trend in a simulation from the five-member ensemble mean, assuming that subannual and interannual to decadal variations in the individual ensemble members are phased randomly and do not imprint on the ensemble mean because they average out. A larger ensemble size would be necessary for this assumption to perfectly hold. However, this potential source of error does not qualitatively alter our results. We remove the seasonal cycle, here defined as the 365 d long mean evolution over the course of a year, from the ensemble mean by smoothing the ensemble mean with a 365 d running mean filter, i.e., by calculating the convolution of the time series with a rectangular window of length 365 and height 1∕365. This filter also removes variability on subannual and interannual timescales and thereby also reduces the error we make due to the small ensemble size that is discussed above. We then subtract the running-mean-filtered ensemble mean from the five ensemble members to remove the secular trend in the individual ensemble members.

Table A1Simulated global ensemble-mean ΩA extreme-event characteristics, when extremes are defined with respect to a shifting baseline. Values in brackets denote ensemble minima and maxima.

Download Print Version | Download XLSX

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f14

Figure A1Simulated globally averaged changes in [H+] extreme events defined with respect to the fixed preindustrial baseline. Shown are changes over the 1861–2100 period following historical (black lines) and future RCP8.5 (red) and RCP2.6 (blue) scenarios for maximal intensity at the surface (a) and at 200 m (b), duration at the surface (c) and at 200 m (d), yearly extreme days at 200 m (e), and volume in the upper 200 m (f). The thick lines display the five-member ensemble means, and the shaded areas represent the maximum and minimum ranges of the individual ensemble members.

Download

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f15

Figure A2Simulated regional changes in [H+] extreme-event characteristics between the preindustrial period and 2081–2100 following the RCP2.6 scenario. The extreme events are defined with respect to shifting baselines. Shown are the changes in yearly extreme days (a, b), maximal intensity (c, d), and duration (e, f). Left panels show changes for the surface, whereas right panels show changes for 200 m. Shown are changes averaged over all five ensemble members. The black contours highlight the pattern structures. Gray areas represent areas with no extremes during 2081–2100.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f16

Figure A3Simulated characteristics of surface [H+]  extreme events for the preindustrial period (a, b), 1986–2005 ensemble mean (c–e), RCP8.5 2081–2100 ensemble mean (f–h), and RCP2.6 2081–2100 ensemble mean (i–k). The extreme events are defined with respect to shifting baselines. Gray colors represent regions where no ensemble member simulates extremes. The black contours highlight the pattern structures.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f17

Figure A4Simulated characteristics of [H+]  extreme events at 200 m for the preindustrial period (a, b), 1986–2005 ensemble mean (c–e), RCP8.5 2081–2100 ensemble mean (f–h), and RCP2.6 2081–2100 ensemble mean (i–k). The extreme events are defined with respect to shifting baselines. Gray colors represent regions where no ensemble member simulates extremes. The black contours highlight the pattern structures.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f18

Figure A5Simulated ensemble-mean changes in mean CT (a, e), AT (b, f), T (c, g), and S (d, h) from the preindustrial period to 2081–2100 following the RCP8.5 scenario. Shown are changes for (a–d) the surface and (e–h) at 200 m. The black contours highlight the pattern structures.

https://bg.copernicus.org/articles/17/4633/2020/bg-17-4633-2020-f19

Figure A6Simulated ensemble-mean changes in the variances of CT (a, e), AT (b, f) T (c, g), and S (d, h) from the preindustrial period to 2081–2100 under the RCP8.5 scenario. Shown are changes for (a–d) the surface and (e–h) at 200 m. The black contours highlight the pattern structures.

Appendix B: Identifying interannual and subannual variability

The spectral density describes how the variance in a time series is distributed over different frequencies νj. It is proportional to the absolute value squared of the discrete Fourier transformation (DFT) of the time series. Defining the spectral density only for positive frequencies, it is given by

(B1) f ( ν j ) = 2 Δ t 2 T k = 1 N x k exp - i 2 π ν j Δ t k 2 ,

with N the number of time steps, xk the values of the time series at each time step, Δt the time interval between two time steps, T=NΔt, and the frequencies νj=j/T. The autocovariance is the inverse Fourier transform of the spectral density (Wiener–Khintchine theorem, Chatfield1996). In the continuous case, the theorem states

(B2) γ ( τ ) = - f ̃ ( ν ) exp ( i 2 π ν τ ) d ν ,

with the autocovariance function γ(τ) and the spectral density f̃ defined for positive and negative frequencies. Since the two-sided spectral density, f̃, is a real and even function, one can also use

(B3) γ ( τ ) = 0 f ( ν ) cos ( 2 π ν τ ) d ν

with the one-sided spectral density f=2f̃ that is used in this text. As a consequence, the variance within the time series, given by the autocovariance at lag zero, is obtained by integrating the spectral density over all positive frequencies, σ2=0f(ν)dν. For a discrete time series, where the maximal resolved frequency is given by νmax=1/2Δt, the identity reads

(B4) σ 2 = j = 0 N / 2 f ( ν j ) 1 N Δ t .

Based on this equation, one can separate the contributions to variance from low-frequency and high-frequency variations. In this study, we determine interannual variability and subannual variability. Interannual variability is calculated by summing over the contributions to variance from all frequencies up to a cycle of once per year, i.e., by evaluating the sum up to icut for which νcut=1/365 d−1. Accordingly, subannual variability is obtained by evaluating the sum from icut+1 to N∕2. Prior to this separation, the seasonal variability is removed from the data by subtracting the 365 d climatology.

Appendix C: Decomposition of [H+] variance change

Following Eq. (2) in the main text, the variance in [H+] (or ΩA) can be approximated as a function of the four sensitivities

(C1) s = H + A T , H + C T , H + S , H + T

that in turn depend on the mean values of the drivers, the four standard deviations of the drivers

(C2) σ = σ A T , σ C T , σ S , σ T ,

and the six pairwise correlation coefficients, in matrix notation given by

(C3) ρ = 1 ρ A C ρ A S ρ A T ρ A C 1 ρ C S ρ C T ρ A S ρ C S 1 ρ S T ρ A T ρ C T ρ S T 1 .

Based on this notation, we can rewrite Eq. (2) of the main text as

(C4) σ H + 2 = i = 1 4 j = 1 4 s i s j σ i σ j ρ i j .

We use Eq. (C4) and decompose the variability change between the preindustrial period and 2081–2100 into the contributions from changes in s, σ, and ρ based on a Taylor expansion. Since [H+] variance represented by Eq. (C4) is a polynomial of fifth order in these variables, its Taylor series has five nonvanishing orders. We use the drivers' standard deviations instead of their variances for the decomposition. With the latter, the Taylor expansion would have infinite terms and could not be decomposed exactly as it is done in the following. However, it would asymptotically lead to the same decomposition of [H+] variance change into ΔsσH+2, ΔσσH+2, ΔsσσH+2, and Δρ+σH+2 that is presented below. Furthermore, it should be noted that the resulting decomposition of [H+] variance change only approximates the simulated variance change because it is based on Eq. (C4), which itself is based on a first-order Taylor expansion of [H+] with respect to the drivers.

In the following, all terms of the Taylor series are given. We denote the sum of first-order terms that contain changes in the four sensitivities Δs1,…4 by Δs(1)σH+2, the sum of second-order terms that contain changes in the sensitivities and standard deviations by Δsσ(2)σH+2, and so on.

The first order is given by Δ(1)σH+2=Δs(1)σH+2+Δσ(1)σH+2+Δρ(1)σH+2 with

(C5) Δ s ( 1 ) σ H + 2 = 2 k = 1 4 j = 1 4 s j σ k σ j ρ k j Δ s k , Δ σ ( 1 ) σ H + 2 = 2 k = 1 4 j = 1 4 s k s j σ j ρ k j Δ σ k , Δ ρ ( 1 ) σ H + 2 = k = 1 4 l = 1 4 s k s l σ k σ l Δ ρ k l .

The second order contains

(C6) Δ s s ( 2 ) σ H + 2 = k = 1 4 l = 1 4 σ k σ l ρ k l Δ s k Δ s l , Δ σ σ ( 2 ) σ H + 2 = k = 1 4 l = 1 4 s k s l ρ k l Δ σ k Δ σ l , Δ s σ ( 2 ) σ H + 2 = 2 k = 1 4 l = 1 4 s l σ l ρ k l Δ s k Δ σ k + s l σ k ρ k l Δ s k Δ σ l , Δ s ρ ( 2 ) σ H + 2 = 2 k = 1 4 l = 1 4 s l σ k σ l Δ s k Δ ρ k l , Δ σ ρ ( 2 ) σ H + 2 = 2 k = 1 4 l = 1 4 s k s l σ l Δ σ k Δ ρ k l .

The third-order terms read

(C7) Δ s s σ ( 3 ) σ H + 2 = 2 k = 1 4 l = 1 4 σ l ρ k l Δ s k Δ s l Δ σ k , Δ s σ σ ( 3 ) σ H + 2 = 2 k = 1 4 l = 1 4 s l ρ k l Δ s k Δ σ k Δ σ l , Δ s s ρ ( 3 ) σ H + 2 = k = 1 4 l = 1 4 σ k σ l Δ s k Δ s l Δ ρ k l , Δ σ σ ρ ( 3 ) σ H + 2 = k = 1 4 l = 1 4 s k s l Δ σ k Δ σ l Δ ρ k l , Δ s σ ρ ( 3 ) σ H + 2 = 2 k = 1 4 l = 1 4 s l σ k Δ s k Δ σ l Δ ρ k l + s l σ l Δ s k Δ σ k Δ ρ k l .

The fourth order reads

(C8) Δ s s σ σ ( 4 ) σ H + 2 = k = 1 4 l = 1 4 ρ k l Δ s k Δ s l Δ σ k Δ σ l , Δ s s σ ρ ( 4 ) σ H + 2 = 2 k = 1 4 l = 1 4 σ l Δ s k Δ s l Δ σ k Δ ρ k l , Δ s σ σ ρ ( 4 ) σ H + 2 = 2 k = 1 4 l = 1 4 s l Δ s k Δ σ k Δ σ l Δ ρ k l .

And the fifth order is given by

(C9) Δ s s σ σ ρ ( 5 ) σ H + 2 = k = 1 4 l = 1 4 Δ s k Δ s l Δ σ k Δ σ l Δ ρ k l .

We identify the variance change from changes in the sensitivities as

(C10) Δ s σ H + 2 = Δ s ( 1 ) σ H + 2 + Δ s s ( 2 ) σ H + 2 ,

the change from standard deviation changes as

(C11) Δ σ σ H + 2 = Δ σ ( 1 ) σ H + 2 + Δ σ σ ( 2 ) σ H + 2 ,

the change from simultaneous changes in sensitivities and standard deviations as

(C12) Δ s σ σ H + 2 = Δ s σ ( 2 ) σ H + 2 + Δ s s σ ( 3 ) σ H + 2 + Δ s σ σ ( 3 ) σ H + 2 + Δ s s σ σ ( 4 ) σ H + 2 ,

and that from correlation changes and mixed contributions that include correlation changes as

(C13) Δ ρ + σ H + 2 = Δ ρ ( 1 ) σ H + 2 + Δ s ρ ( 2 ) σ H + 2 + Δ σ ρ ( 2 ) σ H + 2 + Δ s s ρ ( 3 ) σ H + 2 + Δ σ σ ρ ( 3 ) σ H + 2 + Δ s σ ρ ( 3 ) σ H + 2 + Δ s s σ ρ ( 4 ) σ H + 2 + Δ s σ σ ρ ( 4 ) σ H + 2 + Δ s s σ σ ρ ( 5 ) σ H + 2 .

Finally, we calculate the analogs of Eqs. (C10)–(C13) that only take into account changes in CT; changes in CT and AT; and changes in CT, AT, and T. This is done by calculating Δs1,…4 only based on mean changes in the considered variables and by setting the standard deviation changes in variables and correlation changes in pairs of variables that are not considered to zero.

Appendix D: Comparison of simulated ensemble-mean trends in seasonal amplitude to observation-based trends

We construct confidence intervals for the observation-based slope estimates following Hartmann et al. (2013). For the simulations, we use the arithmetic average of the five ensemble-member slope estimates, b^k, as the estimator,

(D1) b ^ = 1 5 k = 1 5 b ^ k ,

with estimated variance

(D2) σ ^ b 2 = 1 5 2 k = 1 5 σ ^ b k 2 .

We then construct the confidence interval for b^ as

(D3) ( b ^ - q σ ^ b , b ^ + q σ ^ b ) ,

with q the (1+p)/2 quantile (we use p=0.9) of the t distribution with 5(N-2) degrees of freedom. We correct the sample size N (34, the number of years we use for the fits) to a reduced sample size Nr when we find positive lag-one autocorrelation in the residuals of the fits (data – linear regression model). Lag-one autocorrelation is estimated as the average of the five ensemble-member lag-one autocorrelation estimates,

(D4) ρ ^ = 1 5 k = 1 5 ρ ^ k ,

and we obtain Nr=N(ρ^-1)/(ρ^+1). Positive ρ^ is only found in the northern high latitudes. This is in contrast to the observation-based case, where we find large positive ρ^o (up to 0.7) in the residuals of all latitude bands besides the tropical region.

For testing the significance of a difference between the simulation slope estimate b^ and the observation-based estimate b^o, we use Welch's test, which assumes different variances for the two estimates (Andrade and Estévez-Pérez2014). The variance of the simulation slope estimate is calculated by dividing the ensemble-averaged slope variance by the ensemble size (Eq. D2) and is hence smaller than the observation-based slope variance. If the absolute value of the test statistic

(D5) b ^ - b ^ o σ ^ b 2 + σ ^ o

is larger than the (1+p)/2 quantile of the t distribution with (Andrade and Estévez-Pérez2014)

(D6) σ ^ b 2 + σ ^ b o 2 2 σ ^ b 4 / 5 ( N r - 2 ) + σ ^ b o 4 / ( N r , o - 2 )

degrees of freedom, we consider the observation-based and simulation slope to be different from each other with a confidence level of p=0.9.

Data availability

The GFDL ESM2M model data underlying the figures and analyses are available under https://zenodo.org/record/4032577 (Burger2020).

Author contributions

FAB and TLF designed the study. FAB performed the simulations, assisted by TLF and JGJ. FAB performed the analysis and wrote the initial manuscript. All authors contributed to the writing of the paper.

Competing interests

The authors declare that they have no conflict of interest.

Disclaimer

The work reflects only the authors' view; the European Commission and their executive agency are not responsible for any use that may be made of the information the work contains.

Acknowledgements

Friedrich A. Burger and Thomas L. Frölicher thank the CSCS Swiss National Supercomputing Centre for computing resources. The authors thank Elizabeth Drenkard, Fortunat Joos, and Jens Terhaar for discussions and comments; Rick Slater for the help in porting the ESM2M model code to CSCS; and James Orr, Sarah Schlunegger, Jean-Pierre Gattuso, and one anonymous reviewer for their excellent and insightful reviews.

Financial support

This research has been supported by the Swiss National Science Foundation (grant no. PP00P2_170687) and Horizon 2020 (COMFORT, Our common future ocean in the Earth system – quantifying coupled cycles of carbon, oxygen, and nutrients for determining and achieving safe operating spaces with respect to tipping points (grant no. 820989)).

Review statement

This paper was edited by Jean-Pierre Gattuso and reviewed by James Orr, Sarah Schlunegger, and one anonymous referee.

References

Anderson, J. L., Balaji., V., Broccoli, A. J., Cooke, W. F., Delworth, T. L., Dixon, K. W., Donner, L. J., Dunne, K. A., Freidenreich, S. M., Garner, S. T., Gudgel, R. G., Gordon, C. T., Held, I. M., Hemler, R. S., Horowitz, L. W., Klein, S. A., Knutson, T. R., Kushner, P. J., Langenhost, A. R., Cheung, L. N., Liang, Z., Malyshev, S. L., Milly, P. C. D., Nath, M. J., Ploshay, J. J., Ramaswamy, V., Schwarzkopf, M. D., Shevliakova, E., Sirutis, J. J., Soden, B. J., Stern, W. F., Thompson, L. A., Wilson, R. J., Wittenberg, A. T., and Wyman, B. L.: The New GFDL Global Atmosphere and Land Model AM2–LM2: Evaluation with Prescribed SST Simulations, J. Clim., 17, 4641–4673, https://doi.org/10.1175/JCLI-3223.1, 2004. a

Andrade, J. and Estévez-Pérez, M.: Statistical comparison of the slopes of two regression lines: A tutorial, Anal. Chim. Acta, 838, 1–12, https://doi.org/10.1016/j.aca.2014.04.057, 2014. a, b

Bakker, D. C. E., Pfeil, B., O'Brien, K. M., Currie, K. I., Jones, S. D., Landa, C. S., Lauvset, S. K., Metzl, N., Munro, D. R., Nakaoka, S.-I., Olsen, A., Pierrot, D., Saito, S., Smith, K., Sweeney, C., Takahashi, T., Wada, C., Wanninkhof, R., Alin, S. R., Becker, M., Bellerby, R. G. J., Borges, A. V., Boutin, J., Bozec, Y., Burger, E., Cai, W.-J., Castle, R. D., Cosca, C. E., DeGrandpre, M. D., Donnelly, M., Eischeid, G., Feely, R. A., Gkritzalis, T., González-Dávila, M., Goyet, C., Guillot, A., Hardman-Mountford, N. J., Hauck, J., Hoppema, M., Humphreys, M. P., Hunt, C. W., Ibánhez, J. S. P., Ichikawa, T., Ishii, M., Juranek, L. W., Kitidis, V., Körtzinger, A., Koffi, U. K., Kozyr, A., Kuwata, A., Lefèvre, N., Lo Monaco, C., Manke, A., Marrec, P., Mathis, J. T., Millero, F. J., Monacci, N., Monteiro, P. M. S., Murata, A., Newberger, T., Nojiri, Y., Nonaka, I., Omar, A. M., Ono, T., Padín, X. A., Rehder, G., Rutgersson, A., Sabine, C. L., Salisbury, J., Santana-Casiano, J. M., Sasano, D., Schuster, U., Sieger, R., Skjelvan, I., Steinhoff, T., Sullivan, K., Sutherland, S. C., Sutton, A., Tadokoro, K., Telszewski, M., Thomas, H., Tilbrook, B., van Heuven, S., Vandemark, D., Wallace, D. W., and Woosley, R.: Surface Ocean CO2 Atlas (SOCAT) V4, PANGAEA, https://doi.org/10.1594/PANGAEA.866856, 2016. a

Bednaršek, N., Tarling, G. A., Bakker, D. C. E., Fielding, S., Jones, E. M., Venables, H. J., Ward, P., Kuzirian, A., Lézé, B., Feely, R. A., and Murphy, E. J.: Extensive dissolution of live pteropods in the Southern Ocean, Nat. Geosci., 5, 881–885, https://doi.org/10.1038/ngeo1635, 2012. a

Bednaršek, N., Feely, R. A., Reum, J. C. P., Peterson, B., Menkel, J., Alin, S. R., and Hales, B.: Limacina helicina shell dissolution as an indicator of declining habitat suitability owing to ocean acidification in the California Current Ecosystem, P. R. Soc. B, 281, 20140, https://doi.org/10.1098/rspb.2014.0123, 2014. a

Bindoff, N. L., Cheung, W. W. L., Kairo, J. G., Arístegui, J., Guinder, V. A., Hallberg, R., Hilmi, N., Jiao, N., Karim, M. S., Levin, L., O'Donoghue, S., Purca Cuicapusa, S. R., Rinkevich, B., Suga, T., Tagliabue, A., and Williamson, P.: Changing Ocean, Marine Ecosystems, and Dependent Communities, in: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, edited by: Pörtner, H.-O., Roberts, D. C., Masson-Delmotte, V., Zhai, P., Tignor, M., Poloczanska, E., Mintenbeck, K., Alegría, A., Nicolai, M., Okem, A., Petzold, J., Rama, B., and Weyer, N. M., in press, 2020. a, b

Bopp, L., Resplandy, L., Orr, J. C., Doney, S. C., Dunne, J. P., Gehlen, M., Halloran, P., Heinze, C., Ilyina, T., Séférian, R., Tjiputra, J., and Vichi, M.: Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models, Biogeosciences, 10, 6225–6245, https://doi.org/10.5194/bg-10-6225-2013, 2013. a, b

Britton, D., Cornwall, C. E., Revill, A. T., Hurd, C. L., and Johnson, C. R.: Ocean acidification reverses the positive effects of seawater pH fluctuations on growth and photosynthesis of the habitat-forming kelp, Ecklonia radiata, Sci. Rep., 6, 26036, https://doi.org/10.1038/srep26036, 2016. a

Burger, F. A.: Data used for figures and tables in “Increase in ocean acidity variability and extremes under increasing atmospheric CO2”, Zenodo, https://doi.org/10.5281/zenodo.4032577, 2020. a

Caldeira, K. and Wickett, M. E.: Anthropogenic carbon and ocean pH, Nature, 425, 365–365, https://doi.org/10.1038/425365a, 2003. a, b

Carter, B. R., Frölicher, T. L., Dunne, J. P., Rodgers, K. B., Slater, R. D., and Sarmiento, J. L.: When can ocean acidification impacts be detected from decadal alkalinity measurements?, Global Biogeochem. Cy., 30, 595–612, https://doi.org/10.1002/2015GB005308, 2016. a

Carter, B. R., Feely, R. A., Williams, N. L., Dickson, A. G., Fong, M. B., and Takeshita, Y.: Updated methods for global locally interpolated estimation of alkalinity, pH, and nitrate, Limnol. Oceanogr. Method., 16, 119–131, https://doi.org/10.1002/lom3.10232, 2018. a

Chatfield, C.: The analysis of time series – an introduction, Chapman & Hall/CRC, 5th Edn., 1996. a, b

Cheung, W. W. L. and Frölicher, T. L.: Marine heatwaves exacerbate climate change impacts for fisheries in the northeast Pacific, Sci. Rep., 10, 6678, https://doi.org/10.1038/s41598-020-63650-z, 2020. a

Coles, S.: An Introduction to Statistical Modeling of Extreme Values, Springer-Verlag, London, 2001. a, b

Collins, M., Sutherland, M., Bouwer, L., Cheong, S.-M., Frölicher, T., Jacot Des Combes, H., Koll Roxy, M., Losada, I., McInnes, K., Ratter, B., Rivera-Arriaga, E., Susanto, R. D., Swingedouw, D., and Tibig, L.: Extremes, Abrupt Changes and Managing Risk, in: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, edited by: Pörtner, H.-O., Roberts, D. C., Masson-Delmotte, V., Zhai, P., Tignor, M., Poloczanska, E., Mintenbeck, K., Alegría, A., Nicolai, M., Okem, A., Petzold, J., Rama, B., and Weyer, N. M., in press, 2020. a

Cornwall, C. E., Comeau, S., DeCarlo, T. M., Larcombe, E., Moore, B., Giltrow, K., Puerzer, F., D'Alexis, Q., and McCulloch, M. T.: A coralline alga gains tolerance to ocean acidification over multiple generations of exposure, Nat. Clim. Change, 10, 143–146, https://doi.org/10.1038/s41558-019-0681-8, 2020. a

Dickson, A. and Millero, F.: A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media, Deep-Sea Res. Pt. A, 34, 1733–1743, https://doi.org/10.1016/0198-0149(87)90021-5, 1987. a

Dickson, A. and Riley, J.: The effect of analytical error on the evaluation of the components of the aquatic carbon-dioxide system, Mar. Chem., 6, 77–85, https://doi.org/10.1016/0304-4203(78)90008-7, 1978. a

Doney, S. C., Fabry, V. J., Feely, R. A., and Kleypas, J. A.: Ocean Acidification: The Other CO2 Problem, Annu. Rev. Mar. Sci., 1, 169–192, https://doi.org/10.1146/annurev.marine.010908.163834, 2009. a, b

Dunne, J. P., John, J. G., Adcroft, A. J., Griffies, S. M., Hallberg, R. W., Shevliakova, E., Stouffer, R. J., Cooke, W., Dunne, K. A., Harrison, M. J., Krasting, J. P., Malyshev, S. L., Milly, P. C. D., Phillipps, P. J., Sentman, L. T., Samuels, B. L., Spelman, M. J., Winton, M., Wittenberg, A. T., and Zadeh, N.: GFDL’s ESM2 Global Coupled Climate–Carbon Earth System Models, Part I: Physical Formulation and Baseline Simulation Characteristics, J. Clim., 25, 6646–6665, https://doi.org/10.1175/JCLI-D-11-00560.1, 2012. a

Dunne, J. P., John, J. G., Shevliakova, E., Stouffer, R. J., Krasting, J. P., Malyshev, S. L., Milly, P. C. D., Sentman, L. T., Adcroft, A. J., Cooke, W., Dunne, K. A., Griffies, S. M., Hallberg, R. W., Harrison, M. J., Levy, H., Wittenberg, A. T., Phillips, P. J., and Zadeh, N.: GFDL’s ESM2 Global Coupled Climate–Carbon Earth System Models, Part II: Carbon System Formulation and Baseline Simulation Characteristics, J. Clim., 26, 2247–2267, https://doi.org/10.1175/JCLI-D-12-00150.1, 2013. a, b

Fassbender, A. J., Rodgers, K. B., Palevsky, H. I., and Sabine, C. L.: Seasonal Asymmetry in the Evolution of Surface Ocean pCO2 and pH Thermodynamic Drivers and the Influence on Sea-Air CO2 Flux, Global Biogeochem. Cy., 32, 1476–1497, https://doi.org/10.1029/2017GB005855, 2018. a, b

Feely, R., Chris, S., Hernandez-Ayon, J., Ianson, D., and Hales, B.: Evidence for Upwelling of Corrosive “Acidified” Water onto the Continental Shelf, Science, 320, 1490–1492, https://doi.org/10.1126/science.1155676, 2008. a

Fischer, E. M. and Knutti, R.: Anthropogenic contribution to global occurrence of heavy-precipitation and high-temperature extremes, Nat. Clim. Change, 5, 560–564, https://doi.org/10.1038/nclimate2617, 2015. a

Form, A. U. and Riebesell, U.: Acclimation to ocean acidification during long-term CO2 exposure in the cold-water coral Lophelia pertusa, Glob. Change Biol., 18, 843–853, https://doi.org/10.1111/j.1365-2486.2011.02583.x, 2012. a

Friedlingstein, P., Jones, M. W., O'Sullivan, M., Andrew, R. M., Hauck, J., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Le Quéré, C., Bakker, D. C. E., Canadell, J. G., Ciais, P., Jackson, R. B., Anthoni, P., Barbero, L., Bastos, A., Bastrikov, V., Becker, M., Bopp, L., Buitenhuis, E., Chandra, N., Chevallier, F., Chini, L. P., Currie, K. I., Feely, R. A., Gehlen, M., Gilfillan, D., Gkritzalis, T., Goll, D. S., Gruber, N., Gutekunst, S., Harris, I., Haverd, V., Houghton, R. A., Hurtt, G., Ilyina, T., Jain, A. K., Joetzjer, E., Kaplan, J. O., Kato, E., Klein Goldewijk, K., Korsbakken, J. I., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lenton, A., Lienert, S., Lombardozzi, D., Marland, G., McGuire, P. C., Melton, J. R., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S.-I., Neill, C., Omar, A. M., Ono, T., Peregon, A., Pierrot, D., Poulter, B., Rehder, G., Resplandy, L., Robertson, E., Rödenbeck, C., Séférian, R., Schwinger, J., Smith, N., Tans, P. P., Tian, H., Tilbrook, B., Tubiello, F. N., van der Werf, G. R., Wiltshire, A. J., and Zaehle, S.: Global Carbon Budget 2019, Earth Syst. Sci. Data, 11, 1783–1838, https://doi.org/10.5194/essd-11-1783-2019, 2019. a

Frölicher, T. L., Sarmiento, J. L., Paynter, D. J., Dunne, J. P., Krasting, J. P., and Winton, M.: Dominance of the Southern Ocean in Anthropogenic Carbon and Heat Uptake in CMIP5 Models, J. Clim., 28, 862–886, https://doi.org/10.1175/JCLI-D-14-00117.1, 2015. a

Frölicher, T. L., Rodgers, K. B., Stock, C. A., and Cheung, W. W. L.: Sources of uncertainties in 21st century projections of potential ocean ecosystem stressors, Global Biogeochem. Cy., 30, 1224–1243, https://doi.org/10.1002/2015GB005338, 2016. a, b

Frölicher, T. L., Fischer, E. M., and Gruber, N.: Marine heatwaves under global warming, Nature, 560, 360–364, https://doi.org/10.1038/s41586-018-0383-9, 2018. a, b, c, d

Frölicher, T. L., Ramseyer, L., Raible, C. C., Rodgers, K. B., and Dunne, J.: Potential predictability of marine ecosystem drivers, Biogeosciences, 17, 2061–2083, https://doi.org/10.5194/bg-17-2061-2020, 2020. a, b

Gallego, M. A., Timmermann, A., Friedrich, T., and Zeebe, R. E.: Drivers of future seasonal cycle changes in oceanic pCO2, Biogeosciences, 15, 5315–5327, https://doi.org/10.5194/bg-15-5315-2018, 2018. a, b

Gattuso, J.-P. and Buddemeier, R. W.: Calcification and CO2, Nature, 407, 311–313, https://doi.org/10.1038/35030280, 2000. a

Gehlen, M., Gruber, N., Gangstro, R., Bopp, L., and Oschlies, A.: Biogeochemical Consequences of Ocean Acidification and Feedback to the Earth System, chap. 12, in: Ocean Acidification, edited by: Gattuso, J.-P. and Hansson, L., Oxford University Press, 230–248, 2012. a

Good, S. A., Martin, M. J., and Rayner, N. A.: EN4: Quality controlled ocean temperature and salinity profiles and monthly objective analyses with uncertainty estimates, J. Geophys. Res.-Ocean., 118, 6704–6716, https://doi.org/10.1002/2013JC009067, 2013. a

Gray, A. R., Johnson, K. S., Bushinsky, S. M., Riser, S. C., Russell, J. L., Talley, L. D., Wanninkhof, R., Williams, N. L., and Sarmiento, J. L.: Autonomous Biogeochemical Floats Detect Significant Carbon Dioxide Outgassing in the High-Latitude Southern Ocean, Geophys. Res. Lett., 45, 9049–9057, https://doi.org/10.1029/2018GL078013, 2018. a

Griffies, S.: ELEMENTS OF MOM4p1, GFDL ocean group technical report no.6, NOAA/Geophysical Fluid Dynamics Laboratory, Princeton University Forrestal Campus, 201 Forrestal Road, Princeton, NJ 08540-6649, 2009. a

Gruber, N.: Warming up, turning sour, losing breath: ocean biogeochemistry under global change, Philos. T. R. Soc. A, 369, 1980–1996, https://doi.org/10.1098/rsta.2011.0003, 2011. a

Hall-Spencer, J. M., Rodolfo-Metalpa, R., Martin, S., Ransome, E., Fine, M., Turner, S. M., Rowley, S. J., Tedesco, D., and Buia, M.-C.: Volcanic carbon dioxide vents show ecosystem effects of ocean acidification, Nature, 454, 96–99, https://doi.org/10.1038/nature07051, 2008. a

Hartmann, D., Klein Tank, A., Rusticucci, M., Alexander, L., Brönnimann, S., Charabi, Y., Dentener, F., Dlugokencky, E., Easterling, D., Kaplan, A., Soden, B., Thorne, P., Wild, M., and Zhai, P.: Observations: Atmosphere and Surface Supplementary Material, Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, 2013. a

Hauck, J. and Völker, C.: Rising atmospheric CO2 leads to large impact of biology on Southern Ocean CO2 uptake via changes of the Revelle factor, Geophys. Res. Lett., 42, 1459–1464, https://doi.org/10.1002/2015GL063070, 2015. a

Hauri, C., Gruber, N., McDonnell, A. M. P., and Vogt, M.: The intensity, duration, and severity of low aragonite saturation state events on the California continental shelf, Geophys. Res. Lett., 40, 3424–3428, https://doi.org/10.1002/grl.50618, 2013. a, b

Hofmann, G. E., Smith, J. E., Johnson, K. S., Send, U., Levin, L. A., Micheli, F., Paytan, A., Price, N. N., Peterson, B., Takeshita, Y., Matson, P. G., Crook, E. D., Kroeker, K. J., Gambi, M. C., Rivest, E. B., Frieder, C. A., Yu, P. C., and Martz, T. R.: High-Frequency Dynamics of Ocean pH: A Multi-Ecosystem Comparison, PloS One, 6, 1–11, https://doi.org/10.1371/journal.pone.0028983, 2011. a, b, c, d, e, f

Jacox, M. G.: Marine heatwaves in a changing climate, Nature, 571, 485–487, https://doi.org/10.1038/d41586-019-02196-1, 2019. a

Joint, I., Doney, S. C., and Karl, D. M.: Will ocean acidification affect marine microbes?, ISME J., 5, 1–7, https://doi.org/10.1038/ismej.2010.79, 2011. a, b

Jones, C. D., Arora, V., Friedlingstein, P., Bopp, L., Brovkin, V., Dunne, J., Graven, H., Hoffman, F., Ilyina, T., John, J. G., Jung, M., Kawamiya, M., Koven, C., Pongratz, J., Raddatz, T., Randerson, J. T., and Zaehle, S.: C4MIP – The Coupled Climate–Carbon Cycle Model Intercomparison Project: experimental protocol for CMIP6, Geosci. Model Dev., 9, 2853–2880, https://doi.org/10.5194/gmd-9-2853-2016, 2016. a

Keller, K. M., Joos, F., and Raible, C. C.: Time of emergence of trends in ocean biogeochemistry, Biogeosciences, 11, 3647–3659, https://doi.org/10.5194/bg-11-3647-2014, 2014. a

Kroeker, K. J., Micheli, F., Gambi, M. C., and Martz, T. R.: Divergent ecosystem responses within a benthic marine community to ocean acidification, P. Natl. Acad. Sci. USA, 108, 14515–14520, https://doi.org/10.1073/pnas.1107789108, 2011. a

Kroeker, K. J., Kordas, R. L., Crim, R., Hendriks, I. E., Ramajo, L., Singh, G. S., Duarte, C. M., and Gattuso, J.-P.: Impacts of ocean acidification on marine organisms: quantifying sensitivities and interaction with warming, Glob. Change Biol., 19, 1884–1896, https://doi.org/10.1111/gcb.12179, 2013. a

Kroeker, K. J., Bell, L. E., Donham, E. M., Hoshijima, U., Lummis, S., Toy, J. A., and Willis-Norton, E.: Ecological change in dynamic environments: Accounting for temporal environmental variability in studies of ocean change biology, Glob. Change Biol., 26, 54–67, https://doi.org/10.1111/gcb.14868, 2020. a

Kwiatkowski, L. and Orr, J. C.: Diverging seasonal extremes for ocean acidification during the twenty-first century, Nat. Clim. Change, 8, 141–145, https://doi.org/10.1038/s41558-017-0054-0, 2018. a, b, c, d, e, f

Kwiatkowski, L., Gaylord, B., Hill, T., Hosfelt, J., Kroeker, K. J., Nebuchina, Y., Ninokawa, A., Russell, A. D., Rivest, E. B., Sesboüé, M., and Caldeira, K.: Nighttime dissolution in a temperate coastal ocean ecosystem increases under acidification, Sci. Rep., 6, 22984, https://doi.org/10.1038/srep22984, 2016. a

Landschützer, P., Gruber, N., and Bakker, D. C. E.: Decadal variations and trends of the global ocean carbon sink, Global Biogeochem. Cy., 30, 1396–1417, https://doi.org/10.1002/2015GB005359, 2016. a, b

Landschützer, P., Gruber, N., Bakker, D. C. E., Stemmler, I., and Six, K. D.: Strengthening seasonal marine CO2 variations due to increasing atmospheric CO2, Nat. Clim. Change, 8, 146–150, https://doi.org/10.1038/s41558-017-0057-x, 2018. a, b, c, d

Lauvset, S. K., Carter, B. R., Perez, F. F., Jiang, L.-Q., Feely, R. A., Velo, A., and Olsen, A.: Processes Driving Global Interior Ocean pH Distribution, Global Biogeochem. Cy., 34, e2019GB006, https://doi.org/10.1029/2019GB006229, 2020. a

Leinweber, A. and Gruber, N.: Variability and trends of ocean acidification in the Southern California Current System: A time series from Santa Monica Bay, J. Geophys. Res.-Ocean., 118, 3622–3633, https://doi.org/10.1002/jgrc.20259, 2013. a

McNeil, B. I. and Matear, R. J.: Southern Ocean acidification: A tipping point at 450-ppm atmospheric CO2, P. Natl. Acad. Sci. USA, 105, 18860–18864, https://doi.org/10.1073/pnas.0806318105, 2008. a

McNeil, B. I. and Sasse, T. P.: Future ocean hypercapnia driven by anthropogenic amplification of the natural CO2 cycle, Nature, 529, 383–386, https://doi.org/10.1038/nature16156, 2016. a, b

Mehrbach, C., Culberson, C. H., Hawley, J. E., and Pytkowicx, R. M.: MEASUREMENT OF THE APPARENT DISSOCIATION CONSTANTS OF CARBONIC ACID IN SEAWATER AT ATMOSPHERIC PRESSURE 1, Limnol. Oceanogr., 18, 897–907, https://doi.org/10.4319/lo.1973.18.6.0897, 1973. a

Morse, J. and Mackenzie, F.: Geochemistry of Sedimentary Carbonates, Elsevier, Amsterdam, 1990. a

Najjar, R. and Orr, J.: Design of OCMIP-2 simulations of chlorofluorocarbons, the solubility pump and common biogeochemistry, internal ocmip report, LSCE/CEA Saclay, Gif-sur-Yvette, France, 1998. a

Oliver, E. C. J., Donat, M. G., Burrows, M. T., Moore, P. J., Smale, D. A., Alexander, L. V., Benthuysen, J. A., Feng, M., Sen Gupta, A., Hobday, A. J., Holbrook, N. J., Perkins-Kirkpatrick, S. E., Scannell, H. A., Straub, S. C., and Wernberg, T.: Longer and more frequent marine heatwaves over the past century, Nat. Commun., 9, 1324, https://doi.org/10.1038/s41467-018-03732-9, 2018. a

Oliver, E. C. J., Burrows, M. T., Donat, M. G., Sen Gupta, A., Alexander, L. V., Perkins-Kirkpatrick, S. E., Benthuysen, J. A., Hobday, A. J., Holbrook, N. J., Moore, P. J., Thomsen, M. S., Wernberg, T., and Smale, D. A.: Projected Marine Heatwaves in the 21st Century and the Potential for Ecological Impact, Front. Mar. Sci., 6, 1–12, https://doi.org/10.3389/fmars.2019.00734, 2019. a, b, c

Oppenheimer, M., Glavovic, B. C., Hinkel, J., van de Wal, R., Magnan, A. K., Abd-Elgawad, A., Cai, R., CifuentesJara, M., DeConto, R. M., Ghosh, T., Hay, J., Isla, F., Marzeion, B., Meyssignac, B., and Sebesvari, Z.: Sea Level Rise and Implications for Low-Lying Islands, Coasts and Communities, in: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, edited by: Pörtner, H.-O., Roberts, D. C., Masson-Delmotte, V., Zhai, P., Tignor, M., Poloczanska, E., Mintenbeck, K., Alegría, A., Nicolai, M., Okem, A., Petzold, J., Rama, B., and Weyer, N. M., in press, 2020. a

Orr, J. C. and Epitalon, J.-M.: Improved routines to model the ocean carbonate system: mocsy 2.0, Geosci. Model Dev., 8, 485–499, https://doi.org/10.5194/gmd-8-485-2015, 2015. a, b

Orr, J. C., Fabry, V. J., Aumont, O., Bopp, L., Doney, S. C., Feely, R. A., Gnanadesikan, A., Gruber, N., Ishida, A., Joos, F., Key, R. M., Lindsay, K., Maier-Reimer, E., Matear, R., Monfray, P., Mouchet, A., Najjar, R. G., Plattner, G.-K., Rodgers, K. B., Sabine, C. L., Sarmiento, J. L., Schlitzer, R., Slater, R. D., Totterdell, I. J., Weirig, M.-F., Yamanaka, Y., and Yool, A.: Anthropogenic ocean acidification over the twenty-first century and its impact on calcifying organisms, Nature, 437, 681–686, https://doi.org/10.1038/nature04095, 2005. a, b, c

Orr, J. C., Epitalon, J.-M., Dickson, A. G., and Gattuso, J.-P.: Routine uncertainty propagation for the marine carbon dioxide system, Mar. Chem., 207, 84–107, https://doi.org/10.1016/j.marchem.2018.10.006, 2018. a

Palter, J. B., Frölicher, T. L., Paynter, D., and John, J. G.: Climate, ocean circulation, and sea level changes under stabilization and overshoot pathways to 1.5 K warming, Earth Syst. Dynam., 9, 817–828, https://doi.org/10.5194/esd-9-817-2018, 2018. a

Resplandy, L., Séférian, R., and Bopp, L.: Natural variability of CO2 and O2 fluxes: What can we learn from centuries-long climate models simulations?, J. Geophys. Res.-Ocean., 120, 384–404, https://doi.org/10.1002/2014JC010463, 2015. a

Riahi, K., Rao, S., Krey, V., Cho, C., Chirkov, V., Fischer, G., Kindermann, G., Nakicenovic, N., and Rafaj, P.: RCP8.5 – A scenario of comparatively high greenhouse gas emissions, Clim. Change, 109, 33–57, https://doi.org/10.1007/s10584-011-0149-y, 2011. a

Riebesell, U., Zondervan, I., Rost, B., Tortell, P. D., Zeebe, R. E., and Morel, F. M. M.: Reduced calcification of marine plankton in response to increased atmospheric CO2, Nature, 407, 364–367, https://doi.org/10.1038/35030078, 2000. a

Rivest, E. B., Comeau, S., and Cornwall, C. E.: The Role of Natural Variability in Shaping the Response of Coral Reef Organisms to Climate Change, Curr. Clim. Change Rep., 3, 271–281, https://doi.org/10.1007/s40641-017-0082-x, 2017. a

Rodgers, K. B., Sarmiento, J. L., Aumont, O., Crevoisier, C., de Boyer Montégut, C., and Metzl, N.: A wintertime uptake window for anthropogenic CO2 in the North Pacific, Global Biogeochem. Cy., 22, 1–16, https://doi.org/10.1029/2006GB002920, 2008. a

Sarmiento, J. and Gruber, N.: Ocean Biogeochemical Dynamics, Princeton University Press, 2006. a

Seneviratne, S., Nicholls, N., Easterling, D., Goodess, C., Kanae, S., Kossin, J., Luo, Y., Marengo, J., McInnes, K., Rahimi, M., Reichstein, M., Sorteberg, A., Vera, C., and Zhang, X.: Changes in climate extremes and their impacts on the natural physical environment, Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation, A Special Report of Working Groups I and II of the Intergovernmental Panel on Climate Change (IPCC), 2012. a, b

Sentman, L. T., Shevliakova, E., Stouffer, R. J., and Malyshev, S.: Time Scales of Terrestrial Carbon Response Related to Land-Use Application: Implications for Initializing an Earth System Model, Earth Interact., 15, 1–16, https://doi.org/10.1175/2011EI401.1, 2011. a

Shevliakova, E., Pacala, S. W., Malyshev, S., Hurtt, G. C., Milly, P. C. D., Caspersen, J. P., Sentman, L. T., Fisk, J. P., Wirth, C., and Crevoisier, C.: Carbon cycling under 300 years of land use change: Importance of the secondary vegetation sink, Global Biogeochem. Cy., 23, 1–16, https://doi.org/10.1029/2007GB003176, 2009. a

Smale, D. A., Wernberg, T., Oliver, E. C. J., Thomsen, M., Harvey, B. P., Straub, S. C., Burrows, M. T., Alexander, L. V., Benthuysen, J. A., Donat, M. G., Feng, M., Hobday, A. J., Holbrook, N. J., Perkins-Kirkpatrick, S. E., Scannell, H. A., Sen Gupta, A., Payne, B. L., and Moore, P. J.: Marine heatwaves threaten global biodiversity and the provision of ecosystem services, Nat. Clim. Change, 9, 306–312, https://doi.org/10.1038/s41558-019-0412-1, 2019. a

Smith, M. D.: An ecological perspective on extreme climatic events: a synthetic definition and framework to guide future research, J. Ecol., 99, 656–663, https://doi.org/10.1111/j.1365-2745.2011.01798.x, 2011. a

Stephenson, D.: Definition, diagnosis, and origin of extreme weather and climate events, in: Climate Extremes and Society, edited by: Murnane, R. and Diaz, H., Cambridge University Press, 11–23, 2008. a

Terhaar, J., Orr, J. C., Ethé, C., Regnier, P., and Bopp, L.: Simulated Arctic Ocean Response to Doubling of Riverine Carbon and Nutrient Delivery, Global Biogeochem. Cy., 33, 1048–1070, https://doi.org/10.1029/2019GB006200, 2019a. a

Terhaar, J., Orr, J. C., Gehlen, M., Ethé, C., and Bopp, L.: Model constraints on the anthropogenic carbon budget of the Arctic Ocean, Biogeosciences, 16, 2343–2367, https://doi.org/10.5194/bg-16-2343-2019, 2019b. a, b

Terhaar, J., Kwiatkowski, L., and Bopp, L.: Emergent constraint on Arctic Ocean acidification in the twenty-first century, Nature, 582, 379–383, https://doi.org/10.1038/s41586-020-2360-3, 2020. a

Turi, G., Alexander, M., Lovenduski, N. S., Capotondi, A., Scott, J., Stock, C., Dunne, J., John, J., and Jacox, M.: Response of O2 and pH to ENSO in the California Current System in a high-resolution global climate model, Ocean Sci., 14, 69–86, https://doi.org/10.5194/os-14-69-2018, 2018. a

van Heuven, S., Pierrot, D., Rae, J., Lewis, E., and Wallace, D.: MATLAB Program Developed for CO2 System Calculations, available at: http://gts.sourceforge.net/ (last access: 13 February 2019), 2011. a

van Vuuren, D. P., Stehfest, E., den Elzen, M. G. J., Kram, T., van Vliet, J., Deetman, S., Isaac, M., Klein Goldewijk, K., Hof, A., Mendoza Beltran, A., Oostenrijk, R., and van Ruijven, B.: RCP2.6: exploring the possibility to keep global mean temperature increase below 2 C, Climatic Change, 109, 95–116, https://doi.org/10.1007/s10584-011-0152-3, 2011. a

Vogel, M. M., Zscheischler, J., Fischer, E. M., and Seneviratne, S. I.: Development of Future Heatwaves for Different Hazard Thresholds, J. Geophys. Res.-Atmos., 125, e2019JD032, https://doi.org/10.1029/2019JD032070, 2020.  a

Weiss, R.: Carbon dioxide in water and seawater: the solubility of a non-ideal gas, Mar. Chem., 2, 203–215, https://doi.org/10.1016/0304-4203(74)90015-2, 1974. a

Wernberg, T., Bennett, S., Babcock, R. C., de Bettignies, T., Cure, K., Depczynski, M., Dufois, F., Fromont, J., Fulton, C. J., Hovey, R. K., Harvey, E. S., Holmes, T. H., Kendrick, G. A., Radford, B., Santana-Garcon, J., Saunders, B. J., Smale, D. A., Thomsen, M. S., Tuckett, C. A., Tuya, F., Vanderklift, M. A., and Wilson, S.: Climate-driven regime shift of a temperate marine ecosystem, Science, 353, 169–172, https://doi.org/10.1126/science.aad8745, 2016. a

Weyer, N. M. (Ed.): Annex I: Glossary, in: Special Report on Ocean and Cryosphere in a Changing Climate, edited by: Pörtner, H.-O., Roberts, D., Masson-Delmotte, V., and Zhai, P., 677–701, Intergovernmental Panel on Climate Change, Geneva, 2019. a

Winton, M.: A Reformulated Three-Layer Sea Ice Model, J. Atmos. Ocean. Technol., 17, 525–531, https://doi.org/10.1175/1520-0426(2000)017<0525:ARTLSI>2.0.CO;2, 2000. a

Wittenberg, A. T., Rosati, A., Delworth, T. L., Vecchi, G. A., and Zeng, F.: ENSO Modulation: Is It Decadally Predictable?, J. Clim., 27, 2667–2681, https://doi.org/10.1175/JCLI-D-13-00577.1, 2014. a, b

Zscheischler, J. and Seneviratne, S. I.: Dependence of drivers affects risks associated with compound events, Sci. Adv., 3, 1–10, https://doi.org/10.1126/sciadv.1700263, 2017. a

Download
Short summary
Ensemble simulations of an Earth system model reveal that ocean acidity extremes have increased in the past few decades and are projected to increase further in terms of frequency, intensity, duration, and volume extent. The increase is not only caused by the long-term ocean acidification due to the uptake of anthropogenic CO2, but also due to changes in short-term variability. The increase in ocean acidity extremes may enhance the risk of detrimental impacts on marine organisms.
Altmetrics
Final-revised paper
Preprint