Drivers of Surface Ocean Acidity Extremes in an Earth System Model

Oceanic uptake of anthropogenic carbon causes acidification, a process that describes the increase in hydrogen ion concentrations ([H+]) and decrease in calcium carbonate mineral saturation states (Ω). Of particular concern are ocean acidity extreme (OAX) events, which pose a significant threat to many calcifying marine organisms. However, the mechanisms driving such extreme events are not well understood. Here, we use high‐frequency output from a fully coupled Earth system model of all processes that influence the surface ocean temperature and carbon budgets and ultimately [H+] and Ω anomalies to quantify the driving mechanisms of the onset and decline of high [H+] and low Ω extreme events. We show that enhanced temperature plays a crucial role in driving [H+] extremes, with increased net ocean heat uptake being the dominant driver of the event onset in the subtropics. In the mid‐to‐high latitudes, decreased downward vertical diffusion and mixing of warm surface waters during summer, and increased vertical mixing with warm and carbon‐rich subsurface waters during winter are the main drivers of high [H+] extreme event onset. In the tropics, increases in vertical advection of carbon‐rich subsurface waters are the primary driver of the onset of high [H+] extremes. In contrast, low Ω extremes are driven in most regions by increases in surface carbon concentration due to increased vertical mixing with carbon‐rich subsurface waters. Our study highlights the complex interplay between heat and carbon anomalies driving OAX events and provides a first foundation for more accurate prediction of their future evolution.


Introduction
Since the beginning of the industrial era, the ocean has taken up 20 to 30 % of the anthropogenic carbon emissions (Friedlingstein et al., 2022).This uptake has caused changes in ocean chemistry, collectively known as ocean acidification (Caldeira & Wickett, 2003;Doney, Fabry, et al., 2009).Specifically, the pH of the surface ocean has decreased by approximately 0.12 since preindustrial times, corresponding to an increase in hydrogen ion concentration ([H + ]) by 30 % (Jiang et al., 2023).In addition, the concentration of carbonate ions ([CO 2− 3 ]) has decreased by about 16 %, which has resulted in a decrease in the saturation state of calcium carbonate (Ω) (Orr et al., 2005;Jiang et al., 2023).These changes are projected to continue and even accelerate in the future (Orr et al., 2005;Steinacher et al., 2009;Bopp et al., 2013;Kwiatkowski et al., 2020;Canadell et al., 2021).By the end of the 21st century, surface ocean [H + ] is projected to increase by another 4 -150 % and [CO 2− 3 ] concentration is projected to decrease by another 2 -48 %, depending on the future carbon emission scenario (Jiang et al., 2023).These ongoing changes in ocean chemistry are expected to have far-reaching implications for marine organisms and the ser-vices they provide to humanity (Kroeker et al., 2013;Doney et al., 2020;Bindoff et al., 2019).
Extreme variations in ocean acidity, known as OAX events, can amplify the impacts of long-term ocean acidification on marine organisms and ecosystems by pushing them beyond their limits of resilience (e.g., Spisla et al., 2021;Gruber et al., 2021;Bednaršek et al., 2022).These events can cause changes in hydrogen ion concentration ([H + ]) and other carbonate system variables of similar magnitude to those expected from longterm ocean acidification during the 21st century (Hofmann et al., 2011;Leinweber & Gruber, 2013;Desmet et al., 2022), particularly in coastal oceans (Torres et al., 2021).OAX events occur on much shorter timescales and can have detrimental impacts on marine organisms, as demonstrated, for example, by laboratory and field studies that show signs of shell dissolution in calcifying organisms after only a few days in undersaturated calcium carbonate waters (e.g., Bednaršek et al., 2012Bednaršek et al., , 2014)).These findings emphasize the need to consider both short-term and long-term impacts of extreme ocean acidity levels when assessing the health and sustainability of marine ecosystems.
OAX events are projected to become more frequent or even permanent due to longterm ocean acidification by the end of the 21st century (Burger et al., 2020).In addition, short-term departures from normal [H + ] conditions are expected also to become larger in the future, since [H + ] becomes more sensitive to variations in physical and biogeochemical ocean conditions as a consequence of the nonlinear nature of oceanic carbon chemistry (Orr et al., 2018;Fassbender et al., 2018;Kwiatkowski et al., 2023).For example, the frequency of [H + ] extreme events relative to a shifting-mean baseline that includes long-term ocean acidification is projected to increase by a factor of 14 under a high emission scenario by the end of the century (Burger et al., 2020).Such increases in extreme departures may further increase the risk for marine ecosystems under ocean acidification, since ecosystems may be pushed earlier and more frequently beyond their limits of resilience.At the same time, variations in [CO − 3 ] and aragonite saturation state (Ω A ) are expected to become smaller, because [CO − 3 ] and Ω A become less sensitive to variations in physical and biogeochemical ocean conditions (Orr et al., 2018;Burger et al., 2020).
Not only the projections of extreme deviations in [H + ] and Ω A from the long-term mean differ.It is also important to note that these extremes often occur independently from each other.For example, the 2013-2015 marine heatwave in the North Pacific, known as 'the Blob', was associated with extremely high [H + ] conditions (Gruber et al., 2021), but not with extremely low Ω A conditions (Mogen et al., 2022).This difference may be attributed to the fact that distinct drivers can cause [H + ] and Ω A extremes.While high [H + ] levels and low Ω A level may arise from increased dissolved inorganic carbon or decreased alkalinity, high [H + ] levels may also be caused by elevated temperatures (Burger et al., 2022).Furthermore, the drivers determine whether acdification extremes co-occur with extremes in other stressors such as temperature.Understanding when these two types of acidification extremes do not coincide is crucial, particularly if expected impacts are primarily linked to one of the two variables.
Most available studies on OAX events have focused on examining their long-term changes under climate change (Burger et al., 2020;Hauri et al., 2013), as well as on identifying the drivers of the mean seasonal cycle (Hagens & Middelburg, 2016;Xue et al., 2021;Orr et al., 2022) and its changes (Kwiatkowski & Orr, 2018).However, the causes of large deviations in [H + ] or Ω A from their mean seasonal cycles during OAX events are currently unknown.These seasonal anomalies are likely driven by changes in temperature and dissolved inorganic carbon (Deser et al., 2010;Doney, Lima, et al., 2009), which are the most important driving variables.The contribution of different physical and biogeochemical processes, such as air-sea heat and CO 2 exchange and vertical mixing of heat and carbon, to changes in surface heat and carbon and ultimately to extremes in [H + ] or Ω A is currently unknown.A better understanding of these processes is cru-cial for making accurate predictions about the future evolution of OAX events at the regional scale (Burger et al., 2020).
In this study, the drivers of extreme events in [H + ] and Ω A in the global surface ocean are analyzed for the first time.The analysis is based on a pre-industrial control simulation of the GFDL ESM2M Earth system model.It makes use of a suite of model tendency terms for the carbon and temperature budgets that allows to decompose changes in temperature and carbon into contributions from the underlying physical and biogeochemical processes (Gnanadesikan et al., 2012;Palter et al., 2014;S. M. Griffies et al., 2015;Vogt et al., 2022).The remainder of this article is structured as follows.In section 2, the methods used to analyze the drivers of [H + ] extremes are introduced.Section 3 presents the results, and a discussion of the results and conclusions are given in section 4.
We used output of a 100 y preindustrial control simulation that was run under prescribed atmospheric CO 2 levels of 286 ppm (Vogt et al., 2022).Aerosol and solar forcing were also set to preindustrial 1860 values, and no anthropogenic land use and volcanic activity was assumed.We stored output for temperature (T), dissolved inorganic carbon (C T ), total alkalinity (A T ), salinity (S), silicate, and phosphate at two-hourly resolution, which is equivalent to the ocean model time step.By using mocsy 2.0 (Orr & Epitalon, 2015), these data were used to calculate [H + ] and the saturation states of aragonite Ω A -a mineral form of calcium carbonate produced by marine organisms.[H + ] and Ω A were recalculated on the model time step because mocsy 2.0 is also used to calculate partial derivatives of [H + ] and Ω A in the analysis (see section 2.4).This approach thus avoids slight inconsistencies between the carbonate chemistry representations of the ESM2M model and mocsy 2.0, increasing the accuracy of the analysis.The data were then aggregated to daily-mean resolution for the analysis.Additionally, output for the processes that modulate T and C T -specifically T and C T tendency terms -were also stored on two-hourly resolution.Storing tendency terms at each ocean model time step allowed to precisely calculate the changes in daily-mean T and C T arising from the individual tendency terms.

Model evaluation
The findings of this study depend on the models' ability to accurately simulate the variations in [H + ] and Ω A anomalies.The GFDL ESM2M Earth system model, with its nominal 1°horizontal grid resolution in the ocean, is only suitable for assessing ocean acidity extremes on spatial scales of approximately 100 km and larger.The model is not well suited for driver analysis in coastal oceans and at local scales, since mesoscale and submesoscale variability (e.g., Desmet et al., 2022;Hayashida et al., 2020) are not well represented.To evaluate the simulated variability in the open ocean, we compare the model simulation with estimates of observation-based gridded data with a similar 1°horizontal resolution.The observation-based data, covering the period 1982-2021 (see also Burger et al., 2022), consists of the Hadley Centre EN4.2.2 objective analyses T and S fields (Good et al., 2013).Additionally, [H + ] and Ω A were calculated with CO2SYS using SOCAT-based fCO 2 (MPI-SOMFNN v2022; Landschützer et al., 2016;Landschützer et al., 2022) and total alkalinity calculated from S and T using the LIARv2 algorithm (Carter et al., 2018).Since the fCO 2 data is only available on monthly timescales, this model-data comparison is limited to monthly-mean resolution.
After removing the long-term linear trends from the observation-based data, we find a generally good agreement between simulated and observation-based variability of anomalies relative to the seasonal cycle in surface temperature and salinity (Figure 1).
The pattern correlation coefficients of the standard deviation in anomalies are 0.53 for temperature and 0.50 for salinity.However, the model tends to overestimate temperature variability in the Southern Ocean (Figure 1a, e; supporting information Table S3) and salinity variability in the western tropical Pacific and Indian Ocean (Figure 1b, f).S3), particularly in the high latitudes (e.g., +54 % in the Southern Ocean) and the eastern equatorial Pacific.The higher [H + ] variability in the observation-based data is mainly attributable to the historical increase in [H + ] sensitivity with respect to variations in its drivers from ocean acidification (Burger et al., 2020).Recalculating simulated H + ] variability with the driving variables adjusted to the 1982-2021 mean conditions, which include ocean acidification and other historical trends, the excess in observation-based standard deviation is reduced to 4 % globally.However, an excess in observation-based standard deviation of [H + ] anomalies remains in the high latitudes (+18 % over the Southern Ocean).
Calculating observation-based C T following the methodology for [H + ] and Ω A , we find that the remaining mismatch in these regions is associated with a negative bias in simulated variability in C T anomalies (18 % smaller standard deviation in simulated C T seasonal anomalies over the Southern Ocean).It is important to note the uncertainties in the observation-based data from the pCO2 mapping method (Fay et al., 2021), in particular in the high latitudes (Landschützer et al., 2016), highlighting a need to better constrain observation-based carbonate system variability.
In summary, we find a good agreement between simulated and observation-based variability of anomalies relative to the seasonal cycle in all analyzed variables, and a good match in the spatial variability patterns, despite a general low bias in simulated [H + ] variability.These results suggest that the GFDL ESM2M model is well suited to analyze the drivers of extremes in [H + ] and Ω A in the open ocean.

Extreme event definition and identification of onset/decline periods
We examine events of both extremely high [H + ] and extremely low Ω A , which are collectively referred to as OAX events.We define OAX events based on seasonally-varying extreme event thresholds (Hobday et al., 2016;Vogt et al., 2022;Burger et al., 2022).
At each location and for each day of the year, the [H + ] extreme event threshold is determined as the 90th percentile of the 100 anomaly values with respect to the climatological seasonal cycle for that day of the year.As a result, the likelihood that the [H + ] anomaly exceeds the threshold is equal across locations and across the year.The choice of the 90th percentile ensures the inclusion of extreme ocean conditions while maintaining a sufficiently large sample for robust analyses.At a specific location, extreme events in [H + ] are then defined as coherent periods over which the [H + ] anomaly is above the local seasonally varying threshold (Figure 2).Similarly, extremely low Ω A events are defined when Ω A falls below the seasonally-varying thresholds that are given by the 10th percentiles of the anomaly distributions for each calendar day.
At each location and for each OAX event, we identify its onset and decline period (Figure 2).The onset phase is defined as the period between the start of the extreme event (e.g., where the [H + ] anomaly exceeds the seasonally varying threshold) and the peak of the extreme event, where [H + ] anomaly is maximal.Likewise, the decline phase is defined as the period between the peak of the extreme event and the time when [H + ] anomaly falls below the threshold again.In this study, we average the change in [H + ] anomaly and its drivers over these two periods.We assign the day of event peak to the decline period, as the change in [H + ] anomaly on that day characterizes the reduction in [H + ] anomaly between the peak day and the following day.Likewise, the last day of the decline period is excluded, as the change in [H + ] anomaly on that day characterizes the transition from the last day of the event to the first day after the event.

Decomposition of OAX events into drivers
Changes in [H + ] seasonal anomalies (H + ′ ) in each grid cell are decomposed into contributions from T, S, C T and A T (Figure 3; equation ( 1)).The change in [H + ] anomaly between day i and day i+1, denoted by ∆H + ′ (i), is approximated by employing a first order Taylor expansion of [H + ] at day i, and by calculating the seasonal anomalies (denote by primes) of the obtained terms from T, C T , A T , and S: .
(1) ∆T(i), ∆C T (i), ∆A T (i), and ∆S(i) denote the changes in the respective variables between day i and day i + 1.The partial derivatives with respect to T and C T in equation (1) are calculated for each day from daily-mean T, C T , A T , S, silicate and phosphate using mocsy 2.0 (Orr & Epitalon, 2015).The analogous decomposition of anomaly changes is also performed for Ω A .
The approximation of the changes in [H + ] and Ω A seasonal anomalies through the sum of the T, C T , A T , and S terms, as described in equation ( 1) for [H + ], works well.
For example, the root mean squared error (RMSE) over all simulated days of the approximation of [H + ] anomaly change in equation ( 1) is 0.2 pmol kg −1 d −1 (pmol = 10 −12 mol) on global average.RMSE is smaller than 5 % of the standard deviation of [H + ] anomaly change over 99.9 % of the ocean, indicating that the approximation accurately captures variations in [H + ] anomaly change.

Decomposition of T and C T changes during OAX events into tendency terms
Within the ESM2M model, changes in T and C T between two model time steps are calculated from a number of tendencies that describe the changes in T and C T due to the individual physical and biogoechemical processes represented by the model (supporting information text S1; Palter et al., 2014;S. M. Griffies et al., 2015).We make use of these tendency terms to further decompose the changes in T and C T into individual physical and biogeochemical drivers.To do so, changes in daily-mean T or C T due to individual processes are reconstructed by adding up the respective tendency term on two-hourly (model time step) resolution between the two days that are considered (supporting information text S2).1)).The T and CT terms are further decomposed into tendency contributions (equations ( 4) and ( 5)).
For temperature, these individual processes include air-sea exchange of heat (∆T a-s ), resolved and parameterized subgrid-scale horizontal and vertical advection of heat (∆T adv ), vertical diffusion and local mixing of heat (here referred to as vertical diffusion only; ∆T vdiff ), convective vertical mixing of heat in the ocean boundary layer as represented by the nonlocal KPP (K-profile) parametrization (∆T vmix ), and a residual contribution (∆T res ) from other processes, such as neutral diffusion and river runoff (supporting information text S1), as well as grid cell height variations (supporting information text S2): ∆T ≃ ∆T a-s + ∆T vmix + ∆T vdiff + ∆T adv + ∆T res . (2) Likewise, for C T the contributions include air-sea exchange of CO 2 (∆C a-s T ), resolved and parameterized subgrid-scale horizontal and vertical advection of carbon (∆C adv T ), vertical diffusion and local mixing of carbon (∆C vdiff T ), nonlocal KPP convective mixing of carbon (∆C vmix T ), biological carbon uptake and release (∆C bio T ), and other processes including grid cell height variations (∆C res T ): More details on the individual tendencies and their underlying parametrizations can be found in supporting information text S1.The tendencies from grid cell height variations (part of the ∆T res and ∆C res T terms) do not represent physical or biogeochemical processes.However, they are needed to precisely reproduce ∆T and ∆C T with equations (2) and (3).Based on equations ( 2) and ( 3), the T and C T terms in equation ( 1) are decomposed into the individual tendency contributions: . (5) The analogous decomposition is also performed for Ω A .

Results
In Section 3.1, we quantify the contributions of the four drivers T, C T , A T , and S to the onset and decline of surface high [H + ] extremes.Subsequently, we evaluate the specific processes that modulate the impact of the two most important drivers, temperature (Section 3.2) and carbon concentrations (Section 3.3).We investigate the seasonal differences in these processes in Section 3.4.In Section 3.5, we briefly compare our findings on high [H + ] levels with those for low Ω A extremes.(in 94 % of the ocean surface area), with the largest decreases simulated in the tropical Pacific and the high latitudes (Figure 4f).Similarly, the temperature term also decreases in most regions (in 92 % of the surface ocean), with the most pronounced decreases in the subtropics (Figure 4d), where the temperature term is the main driver of [H + ] decline.An exception is again the equatorial Pacific, where temperature increases during the decline period of [H + ] extremes, thereby counteracting event decline.

Drivers of temperature variations during [H + ] extreme events
To understand the individual processes driving the changes in temperature anomalies during [H + ] extreme events and hence the temperature contribution to onset and decline of [H + ] extreme events (Figure 4c, d), the change in temperature anomaly during event onset and decline is decomposed (see equation ( 4); Table 1) into the contributions from air-sea heat exchange (Figure 5b, c), nonlocal KPP convective mixing (Figure 5e, f), vertical diffusion and local mixing (Figure 5h, i), and horizontal and vertical advection (Figure 5k, l).For ease of interpretation of the anomaly patterns, we also show the climatological means of the tendency contributions to the T term (Figure 5a,d,g,j).These are determined by calculating the temporal mean values of the tendency contributions to the T term instead of calculating their seasonal anomalies.
At the global scale, reduced ocean heat loss (i.e., net ocean heat uptake) contributes most to the increases in temperature anomalies during the onset period of [H + ] extreme events (Table 1).The net ocean heat uptake increases At the regional scale (Figure 5), the positive contribution from air-sea heat exchange is largest in the low-to-mid latitudes and in particular in the subtropical oceans (Table 1; Figure 5b), while the contribution is much smaller or negative in the high latitudes.In the subtropics, air-sea heat exchange often changes the sign from net loss to the atmosphere to net uptake during [H + ] extreme events (Figure 5a,b).Vertical diffusion and local mixing increases temperature anomaly and hence contributes positively to the T term in all ocean regions, except in the tropical Pacific and Indian Ocean, where vertical diffusion and local mixing of heat to the subsurface is increased during [H + ] extreme events (Figure 5h).The vertical diffusion and local mixing contribution is most positive in the North Atlantic, North Pacific, and Southern Ocean (Table 1, Figure 5h).In the subtropics and the mid-to-high latitudes during summer, the positive contribution from vertical diffusion and local mixing arises due to a reduction in mixing with colder subsurface waters that coincides with a reduction in wind strength (supporting information Figure S3a, c).In the high latitudes during winter, the positive contribution instead often arises due to an increase in upward mixing of heat.The increase in upward mixing of heat is associated with stronger winds (supporting information Figure S3c).Decreases in nonlocal KPP convective mixing decrease surface temperature and [H + ] almost in the entire global surface ocean, especially in the subtropics (Table 1, Figure 5e), where also increases in temperature anomaly due to air-sea heat exchange are largest (Figure 2b, Figure 5b).The contribution from advective heat transport is generally small (Figure 5k), except in the tropical Pacific, where it decreases temperature and [H + ] during the onset of [H + ] extreme events (Table 1, Figure 5e).During the decline phase of [H + ] extreme events, the temperature decrease in the subtropics mainly results from increased heat losses to the atmosphere (Figure 5c).Vertical diffusion and local mixing also decreases temperature over most of the ocean (Fig- ure 5i) and is the main driver of temperature decrease in the Southern Ocean (Table 1).
The increases in temperature that counteract the [H + ] event decline in the tropical Pacific (Figure 4d) result from enhanced ocean heat uptake during the decline phase (Fig- ure 5c).

Drivers of carbon variations during [H + ] extreme events
At the global scale ( ology are small at the global scale.The residual term for C T is larger than for T and mainly stems from neutral diffusion and a tendency that compensates a numerical artifact associated with the smoothing of the free ocean surface in the model.At the regional scale (Figure 6), the contribution from vertical and horizontal advection to the onset of [H + ] extreme events is largest in the tropics (Figure 6k, Table 1).
Smaller positive contributions from advection are also simulated in high-latitude regions.
However, the advection contribution is slightly negative in the subtropics.In the ESM2M  creasing C T ) in the low-to-mid latitudes coincide with negative anomalies in wind stress during event onset (supporting information Figure S3a, c).In the high-latitude regions where positive anomalies in vertical diffusion and local mixing of temperature and C T are simulated (Figures 5h and 6h), also wind stress and mixed layer depth (not shown) are increased during event onset, in particular during winter (supporting information Fig- ure S3c).The increased wind stress may be the reason for enhanced vertical mixing during event onset in these regions.
Biological activity generally reduces C T everywhere, because biological production outweighs decomposition at the surface (Figure 6m).During the onset of [H + ] extreme events, increases in biological production decrease [H + ] anomaly in the tropics (Table 1), while reductions in biological production increase [H + ] anomaly in the mid-to-high latitudes (Figure 6n).In the tropical regions, increased nutrient concentrations are simulated during OAX events (not shown), which may cause increased phytoplankton growth there.In the mid-to-high latitudes, low biological production may be connected to nutrient limitation and / or low phytoplankton biomass due to enhanced zooplankton grazing under the elevated temperatures.
During the decline phase of [H + ] extremes, C T anomaly decreases almost in the entire ocean (Figure 4f).This decrease is mainly due to loss of carbon to the atmosphere (Figure 6c), which remains similarly strong as during the onset phase.In the tropical ocean, biological production continues to decrease C T anomaly during event decline (Fig- ure 6o), and also advection reduces C T anomaly during event decline there (Figure 6c).
The convective mixing term balances the carbon losses from air-sea gas exchange and counteracts [H + ] event decline everywhere in the ocean (Figure 6f), and also vertical diffusion and local mixing increases [H + ] anomaly during event decline in most regions (Figure 6i).These increases in vertical mixing, simultaneously also causing temperature decreases in the low-to-mid latitudes (Figure 5i), may be connected to the anomalous heat losses to the atmosphere (Figure 5c), causing loss of buoyancy.

Seasonal variations in drivers of [H + ] extreme event onset
Next, we analyze if the drivers of [H + ] extreme event onset differ between summer and winter season.At the global scale, the changes in [H + ] anomaly per day during extreme event onset are by a factor of two larger during hemispheric summer (April to September on the northern and October to March on the southern hemisphere) than during hemispheric winter (October to March on the northern and April to September on the southern hemisphere; grey bars in Figure 7a).This difference arises from distinct drivers of [H + ] extreme events during individual seasons.Globally, temperature changes are the dominant driver for event onset during the summer months, with 11.8 pmol kg −1 d −1 for temperature vs 1.0 pmol kg −1 d −1 for C T .In contrast, C T changes become more important in winter, with 3.9 pmol kg −1 d −1 for C T vs 2.7 pmol kg −1 d −1 for temperature (Fig- ure 7a).
In the tropics, temperature and C T are both important drivers during hemispheric summer, but C T is the dominant driver during hemispheric winter (Figure 7b).The smaller contribution from temperature in winter is due to lower net ocean heat uptake and increased heat loss from vertical diffusion and local mixing and advection.The larger contribution from C T in winter, on the other hand, is due to larger surface C T increases from advection as well as vertical diffusion and local mixing (supporting information Figure S4b).
In the subtropics, temperature is the dominant driver of [H + ] event onset throughout the year, while the contributions from the other drivers are negligible both in summer and winter (Figure 7c).The overall larger [H + ] increases during event onset in summer result from a larger reduction in heat loss from vertical diffusion and local mixing in summer (supporting information Figure S4c).Results are shown separately for hemispheric summer (April to September on the northern hemisphere and October to March on the southern hemisphere; light colors) and hemispheric winter (remaining months; dark colors).The definition of regions is shown in Figure 4a.Errors of the decomposition are not shown as they are very small (see Table 1).
In the Southern Ocean, the drivers during summer show similarities to subtropical regions, characterized by a large positive temperature contribution and a much smaller and negative C T contribution (Figure 7d).The large positive temperature contribution is due to reduced vertical diffusion and local mixing with colder subsurface waters, associated with reduced winds (supporting information Figure S3a), and due to net ocean heat uptake.The negative C T contribution is caused by reduced mixing with carbonrich subsurface layers and due to air-sea CO 2 loss (supporting information Figure S4d).
The regime is distinctly different during winter, where carbon increases are the main driver of event onset while temperature increases are only of secondary importance (Figure 7d).
The increase in carbon is mainly caused by an increase in vertical diffusion and local mixing with carbon-rich subsurface layers that is associated to enhanced winds (supporting information Figure S3c).It is counteracted by amplified air-sea CO 2 loss to the atmosphere.The smaller positive contribution from temperature increases in winter is caused by enhanced vertical mixing and diffusion, transporting heat from the warmer subsurface to the surface (see also Section 3.2; supporting information Figure S4d).

Drivers of extremes in Ω A
At the global scale, the onset of surface low Ω A extremes is mainly caused by an increase in C T , which accounts for 65 % of the total decrease in Ω A (Figure 8a,b; supporting information Table S1).This is in contrast to high [H + ] extremes, for which temperature accounts for 80 % of the total increase in [H + ].The increase in C T during low Ω A extremes is primarily due to increased local vertical mixing and diffusion (Figure 8d), which is counterbalanced by enhanced biological activity, which decreases C T (Figure 8f), as well as anomalous outgassing of CO 2 (Figure 8c).In addition, the decrease in total alkalinity (Figure 8g) and decrease in temperature from air-sea heat loss (Figure 8h,i) contribute to the decrease in Ω A during event onset.At the regional scale, the primary driver of the onset of low Ω A extremes is the increase in C T over most of the tropics and the Southern Ocean (Figure 8b; supporting information Table S1).The increase in C T contributes 71% to the total decrease in Ω A in the tropics and 99 % in the Southern Ocean.These increases in C T are caused by enhanced vertical mixing and diffusion (Figure 8d), as well as advection of carbon-rich waters (Figure 8e).However, these increases in C T in the tropics and the Southern Ocean are somewhat counterbalanced by decreases in carbon from anomalous CO 2 outgassing (Figure 8c).Enhanced biological activity also decreases carbon in the tropics (Figure 8f).
In the western tropical Pacific, the decrease in C T is driven by the decrease in A T (Fig- ure 8g) resulting from enhanced precipitation (supporting information Figure S9b).Decreases in A T are also contributing to the decreases in Ω A in the subtropical regions, presumably due to enhanced vertical mixing with low-A T subsurface waters.In addition, decreases in temperature from air-sea heat loss contribute to the decreases in Ω A in the subtropics (Figure 8i).
The global decline of surface low Ω A is mainly caused by decreases in C T (predominantly from biological uptake of carbon), increases in A T , and a smaller contribution from increasing temperatures that are caused by surface warming (supporting information Table S1).

Discussion and conclusions
We provide a first assessment of the drivers of surface high [H + ] and low Ω A extreme events using high-frequency output of a comprehensive Earth system model.The results of this modeling study suggest that rising temperatures from enhanced net ocean heat uptake and reduced heat loss through vertical mixing are the primary drivers of high [H + ] events in the subtropical regions and mid-to-high latitudes during summer.In tropical regions, simulated high [H + ] events as well as low Ω A events are often driven by increases in dissolved inorganic carbon due to advection.In mid-to-high latitudes during winter, we also find increased vertical mixing with carbon-rich and often warmer subsurface waters to be an important factor for the onset of high [H + ] events.
Recent studies have investigated the biogeochemical imprint of the 2013-2015 marine heatwave in the North Pacific, known as the Blob.While extremely high [H + ] was identified, the levels of Ω A were not lower than usual, despite both being referred to as OAX events (Gruber et al., 2021;Mogen et al., 2022).Our study offers an explanation for this apparent contradiction: high [H + ] events are often driven by temperature increases, in particular in the North Pacific region where the Blob occurred (Figure 3c), which suggests that the Blob was indeed an [H + ] extreme.On the other hand, low Ω A events usually coincide with periods of decreasing temperature (Figure 7e) and thus unlikely to coincide with MHWs.
The main driver of low Ω A extremes in the model is an increase in C T resulting from vertical mixing, diffusion, and advection, particularly in the tropics and the midto-high latitudes (Figure 7d,e).Similarly, advective increases in C T were also identified as the main driver of [H + ] in the tropics, and vertical diffusion and local mixing was identified as an important driver of [H + ] extremes in the mid-to-high latitudes, in particular during winter (section 3.3).The resemblance in the driving mechanisms in these regions reflects in relatively frequent co-occurrence of simulated high [H + ] and low Ω A extremes: 26 % of high [H + ] extreme days overlap with low Ω A extreme days in the tropics and 31 % of event days overlap in the Southern Ocean (Figure 9), in particular during winter where 43 % of event days coincide.In contrast, the driving mechanisms diverge in the subtropics, reflected in a relatively low overlap of event days of only 11 % there (Figure 9).
Our study suggests that temperature increase is the main driver of [H + ] extremes in the subtropics, which may share similar physical drivers to those of marine heatwaves.
Using the same GFDL ESM2M model, Vogt et al. (2022) found that air-sea heat fluxes were the main factor responsible for temperature increases during the onset of marine heatwaves, particularly in the subtropical oceans, offset by temperature decreases resulting from reduced nonlocal-KPP convective mixing (Figure 1 in Vogt et al., 2022).Our findings are thus consistent with those of Vogt et al. (2022) for the subtropical oceans.
During the onset period, anomalous air-sea heat flux was identified as the primary driver of [H + ] increase, and likewise, during event decline, air-sea heat flux was the main driver of [H + ] decrease (Figure 5b, c; Table 1).Moreover, reduced convective mixing from the nonlocal KPP parameterization was identified as the primary inhibiting factor during the onset period, and to a lesser extent, also an important factor for event decline (Figure 5k, l; Table 1).Given the similarity in the main drivers of [H + ] extremes and marine heatwaves in the subtropical ocean, one can expect that these two univariate extreme events often co-occur.This is indeed the case.S10).Therefore, this study's results on the driving mechanisms of [H + ] extreme event onset often also apply to the preceding preconditioning phase.
Even though we consider our results as robust, a number of caveats need to be discussed.The analysis of OAX event drivers relies on data from an Earth system model, as certain processes cannot be independently validated with observational-based data due to its limited availability.The robustness of our results depends therefore on the Earth system model's accuracy in simulating the physical and biogeochemical processes that precipitation may be less important drivers of [H + ] and Ω A extremes in these regions than identified here.To better constrain the simulated physical and biogeochemical processes, it would be beneficial to compare the identified drivers for the GFDL ESM2M model to those from other Earth system models that can provide the required diagnostic output.
Another shortcoming is that the ocean model used in this study has a relatively coarse spatial resolution and cannot explicitly simulate small-scale circulation features, such as meso-and submesoscale dynamics (S.M. Griffies et al., 2015).Our analysis may therefore underestimate the impact of these small-scale circulation features on the onset and decline of OAX events.Additionally, the coarse resolution of the ocean model limits our analysis to the open ocean, and higher resolution ocean model including improved biogeochemistry would be needed to more accurately represent the drivers in coastal areas (Turi et al., 2018;Terhaar et al., 2019).It should be also noted that the present study focused on the analysis of the mean driving processes over the onset and decline periods of OAX events.However, individual extreme events may be governed by different processes.For example, the drivers in a region can vary between seasons (discussed in Section 3.4), and different types of extremes, characterized by different combinations of drivers, can also occur during the same season (Vogt et al., 2022).The mean process contributions to OAX event onset and decline shown in this study therefore characterize average extremes event in a region and season.Finally, this study analyzes the drivers of OAX events under preindustrial stationary climate conditions.However, ongoing ocean warming and acidification may modify the primary drivers of OAX events, as the background ocean carbon and temperature fields on which the drivers act, as well as the drivers themselves, may change with climate change.To address this limitation, future research should extend our analysis to simulations that include the climate change signal.
In conclusion, our modeling results highlight the crucial role of temperature in driving [H + ] extremes, particularly during summer and in the subtropical oceans.This is primarily attributed to anomalous air-sea heat fluxes and vertical mixing of heat.Furthermore, our results indicate that changes in dissolved inorganic carbon are the dominant driver of low Ω A extremes, as well as for high [H + ] extremes in equatorial regions and high latitudes during winter, thereby designating these regions as hotspots of compound high [H + ] and low Ω A extremes.Our findings enhance our current understanding of OAX events in the global ocean and provide a first foundation for regional predictions of such events.

Open Research Section
The GFDL ESM2M model output and analysis scripts used in this study are available under the Zenodo repository (URL), Friedrich A. Burger (DATE OF PUBLICA-TION).(the repository is underway and the URL and date of publication will be added at the next stage of the review process) module.The Modular Ocean Model version 4p1 (MOM4p1) uses a grid with a horizontal nominal 1 • resolution that increases near the equator to 0.3 • and with a time-varying vertical resolution of about 10 m in the upper ocean.In this study, we analyze data for the uppermost vertical layer that extends from the surface to about 10 m depth.MOM4p1 is coupled to the ocean biogeochemistry model Tracers of Ocean Phytoplankton with Allometric Zooplankton version two (TOPAZv2; Dunne et al., 2013).

Figure 1 .
Figure 1.Standard deviation for anomalies relative to the seasonal cycle in (a,e) T, (b,f) S, (c,g) [H + ], and (d,h) ΩA, of the pre-industrial GFDL ESM2M model simulation (top) and observation-based data over the period 1982-2021 (bottom).The observation-based data was linearly detrended prior to the analysis.

Figure 2 .
Figure 2.An exemplary [H + ] extreme event in the northern subpolar Pacific depicting the event definition as well as the separation into event onset and decline periods.

Figure 3 .
Figure 3.A scheme depicting the decomposition of [H + ] anomaly change (∆H + ′ ) into the T, CT, AT, and S terms (equation (1)).The T and CT terms are further decomposed into tendency

3. 1
Contributions of temperature, carbon, alkalinity and salinity anomalies to the onset and decline of surface [H + ] extremes On a global scale, the average increase in [H + ] anomalies during event onset is 9.1 pmol kg −1 d −1 (Table 1).The main factor contributing to the increase is the rise in temperature during the onset of the event.On a global scale, the temperature increase contributes 7.3 pmol kg −1 d −1 or 80% to the total increase.Increased temperature directly leads to an increase in [H + ] via changes in the carbonate chemistry equilibrium(Zeebe & Wolf-Gladrow, 2001).Increases in C T also contribute to the increase in [H + ] globally, but the contribution of 2.4 pmol kg −1 d −1 is relatively small and accounts for 27% of the total [H + ] increase.Increases in alkalinity slightly counteract the [H + ] increases (-8 %) and contributions from changes in salinity are minor (1 %).At the regional scale (Figure4), we find that the increase in temperature is the dominant driver of [H + ] increases in 78 % of the global ocean surface area during the onset phase, whereas C T dominates over 22 % of the ocean surface area.Regions that are particularly dominated by the C T contribution are the eastern and central tropical Pacific and the Arctic Ocean.There, increases in [H + ] during the onset period result from increases in C T (Figure4e), while temperature decreases make the temperature contribution negative (Figure4c).In the subtropics and the Southern Ocean, increases in temperature and associated increases in positive [H + ] anomalies are somewhat damped by decreases in C T and associated decreases in [H + ] anomalies.In the Kuroshio and Gulf Stream regions and the Northern Indian Ocean, increases in [H + ] anomalies result from increases in both T (Figure4c) and in C T (Figure4e).This pattern also holds true near Antarctica.However, it is important to note that the model may not adequately simulate physical and biogeochemical dynamics close to the Antarctic continent.The contribution from A T is generally smaller and predominately opposite to the contribution of C T (pattern correlation coefficient of -0.73).The A T contribution is most important in the tropical regions where increases in A T tend to decrease [H + ] and therefore inhibit event onset.During the decline period of [H + ] extreme events, simulated reductions in temperature (Figure4d) and C T concentrations (Figure4f) contribute similarly to the decline of [H + ] at the global scale, with a decrease of approximately -4.7 and -5.0 pmol kg −1 d −1 , respectively.This is in contrast to the onset period, where the temperature term dominates at the global scale.At the regional scale, the C T term decreases almost everywhere in pmol kg −1 d −1 ] anomalies during the onset and decline periods of [H + ] extreme events and the tendency contributions to these changes.Results are shown for the global ocean, the tropics, the subtropics, and the Southern Ocean.The regions are defined in Figure4a.For each period and region, the two largest positive and negative contributions are highlighted in red and blue bold font, respectively.The T and CT terms are decomposed into air-sea flux (a-s), convective vertical mixing (vmix ), local vertical mixing and diffusion (vdiff ), advection (adv ), biology (bio), and residual (res) terms.'Total' denotes the sum of all tendency contributions and 'Simulated' denotes the actual simulated [H + ] change during the onset and decline phases.

Figure 4 .
Figure 4.The simulated [H + ] anomaly change during onset and decline of extreme high [H + ] events and the contributions from the T, CT, and AT terms.(a, b) The simulated change in [H + ] anomalies during the onset and decline phases of extreme [H + ] events, the contribution of the (c, d) T term, (e, f) CT term, and (g, h) AT term (see equation (1)).The salinity term is small and not shown (see also Table 1).The solid, dashed, and dotted boxes in a) indicate the tropics (10 • S -10 • N and 220 • W -85 • W in the Pacific, 55 • W -10 • E in the Atlantic, and 50 • E -100 • E in the Indian Ocean), the subtropics (15 • N -30 • N and 205 • W -125 • W in the North Pacific, 30 • S -15 • S and 190 • W -90 • W in the South Pacific, 15 • N -30 • N and 65 • W -25 • W in the North Atlantic, 30 • S -15 • S and 35 • W -5 • E in the South Atlantic, and 30 • S -15 • S and 55 • E -105 • E in the Indian Ocean), and the Southern Ocean (65 • S -45 • S), respectively.These regions are used in Table1.
[H + ] anomalies by 12.6 pmol kg −1 d −1 (138 % of [H + ] anomaly increase) at the global scale.In addition, increases in temperature anomaly associated with reduced vertical diffusion and local mixing of warm waters to the subsurface cause an increase in [H + ] anomalies of 5.0 pmol kg −1 d −1 (54¸,% of [H + ] anomaly increase) during the onset period.Convective mixing increases sea surface temperature by transporting heat to the surface when surface waters lose buoyancy due to heat loss to the atmosphere.This mechanism is less active during the positive airsea heat flux anomalies in the onset period.The associated negative anomalies in nonlocal KPP convective mixing during the onset period reduce surface temperature anomalies and therefore [H + ] anomalies during the onset, strongly dampening the temperatureinduced increases in [H + ] anomalies (−9.5 pmol kg −1 d −1 , -104 % of [H + ] anomaly increase).

Figure 5 .
Figure 5.The decomposition of the T term into tendency contributions.The climatological means of the tendency contributions to the T term (first column) as well as their contributions to the onset (second column) and decline (third column) means of the T term (Figure 4 c, d).
model, C T changes due to advection also include the diluting or concentrating effect on C T from precipitation minus evaporation (Supporting information text S3).Negative precipitation minus evaporation anomalies (i.e., more evaporation than precipitation) increase C T during event onset in the western tropical Pacific and Indian Ocean (Figure S4e), while the advective increases in C T in the remaining tropics result from oceanic advection such as upwelling (supporting information Figure S4h), in particular in the eastern tropical Pacific where also advective decreases in temperature are simulated (Figure 5k).The slightly negative contribution from advection in the subtropics is caused by positive precipitation minus evaporation anomalies during the onset of [H + ] extreme events.The negative anomalies in air-sea CO 2 flux are largest in the high latitudes (Figure 6b).The anomalies in air-sea CO 2 exchange are offset by opposing tendencies from nonlocal KPP convective mixing of carbon in most regions (Figure 6e).The convective mixing increases C T and [H + ] everywhere except in the western tropical Pacific and the tropical Indian Ocean.Vertical diffusion and local mixing generally increase surface C T in the climatological mean (Figure 6g).During the onset of [H + ] extreme events, negative anomalies in vertical diffusion and local mixing counteract increases in [H + ] anomaly in the subtropics, the mid latitudes, and most tropical regions (Figure 6h).In contrast, vertical diffusion and local mixing increases [H + ] anomalies in the eastern tropical Pacific, and in high-latitude regions of the the North Pacific, North Atlantic, and Southern Ocean.Its contribution tends to be opposite to that of temperature vertical diffusion and local mixing due to opposite vertical gradients in temperature and C T .This is not the case in the high-latitude regions where temperature and C T vertical gradients are often both positive towards depth during the winter months.The reductions in vertical diffusion and local mixing of temperature and C T (increasing temperature and de-

Figure 6 .
Figure 6.The decomposition of the CT term into tendency contributions.The climatological means of the tendency contributions to the CT term (first column) as well as their contributions to the onset (second column) and decline (third column) means of the CT term (Figure 4 e, f).

Figure 7 .
Figure 7.The contributions to [H + ] extreme event onset (light and dark grey bars) from T (light and dark orange), CT (light and dark green), and AT and S (light and dark purple).

Figure 8 .
Figure 8.The ΩA anomaly change during the onset phase of low ΩA events, the contributions from the T, CT, and AT terms, and the most important tendency contributions to these.(a) The simulated change in ΩA anomalies during the onset phase.(b-f) The contribution from the CT term (b) and its air-sea CO2 exchange (c), local vertical mixing and diffusion (d), advection (e), and biology (f) contributions.(g) The AT term, and (h) the T term and its contribution from air-sea heat exchange (i).The salinity term and the other tendency contributions are smaller and not shown (see also supporting information TableS1).Blue colors indicate a decrease in ΩA and thus an intensification of low ΩA extremes.

Figure 9 .
Figure 9. Overlap percentage of high [H + ] and low ΩA extremes.It is calculated as the number of days with co-occurring high [H + ] and low ΩA extremes divided by the number of days with high [H + ] or low ΩA extremes (10 % of days).Cold and warm colors indicate that the events coincide less and more frequent than by chance, respectively.
lead to [H + ] and Ω A variations and extremes.As these processes also drive spatial variability patterns in [H + ] and Ω A , a first step is to evaluate simulated variability patterns in the seasonal anomalies of [H + ], Ω A , and the underlying physical fields T and S against observation-based data (Section 2.2).Overall, we found a good agreement both in the magnitude of variability and in spatial differences for all variables.However, simulated variability in [H + ] anomalies in high latitude regions is lower than in the observationbased data, associated with a low bias in simulated variability of C T anomalies.The identified lack in simulated C T variability in these regions suggests a too small contribution from C T variations and a too large contribution from temperature variations to [H + ] dynamics and thus onset and decline of [H + ] extremes (also found for CMIP6-type models inBurger et al. 2022), regionally reinforced by a positive bias in simulated temperature variability in the Southern Ocean.The identified drivers of Ω A extremes are less affected by these biases, since Ω A is less dependent on the balance between temperature and C T variability.Furthermore, simulated salinity and thus freshwater variations are too large in the western tropical Pacific and Indian Ocean.As a result, evaporation and

Table 1 .
The simulated average daily changes in [H +

Table 1 )
, vertical and horizontal advection is the most important driver of C T increase during the onset of [H + ] extreme events, increasing [H + ] anomaly by 5.0 pmol kg −1 d −1 (55 % of [H + ] anomaly increase).In addition, reduced nonlocal KPP convective vertical mixing of carbon increases [H + ] by 2.3 pmol kg −1 d −1 (25 % of [H + ] anomaly increase).These increases are balanced by decreases in C T anomalies from negative anomalies in air-sea CO 2 flux during the onset of [H + ] extreme events (-4.2 pmol kg −1 d −1 ,-46 % of [H + ] anomaly increase).Negative anomalies in air-sea CO 2 flux, i.e., increased carbon loss to the atmosphere or decreased CO 2 uptake from the atmosphere (Figure6a), occur when partial pressure of CO 2 (pCO 2 ) in the surface water is increased.Due to the high correlation between [H + ] and pCO 2 anomalies (Pearson correlation coefficient of 0.99 on global average in the model), negative anomalies in air-sea CO 2 flux during high [H + ] events are expected.The contributions from vertical diffusion, local mixing and bi- Burger et al. (2022)found observationbased evidence for this co-occurrence, indicating that the subtropical oceans are a hotspot for compound high [H + ] and high temperature extremes.The processes responsible for preconditioning [H + ] extremes, i.e., for increasing [H + ] anomaly before crossing the 90th percentile threshold, have not been analyzed yet.However, the simulated drivers during this preconditioning phase are often similar to those during event onset, with temperature increases from air-sea heat flux dominating in the subtropical oceans and advective increases in carbon dominating in the tropical oceans (supporting information TableS2and supporting information Figure