Models project that with current CO 2 emission rates, the Southern Ocean surface will be undersaturated with respect to aragonite by the end of this century [1][2][3][4] . This will result in widespread impacts on biogeochemistry and ocean ecosystems [5][6][7] , particularly the health of aragonitic organisms, such as pteropods 7 , which can dominate polar surface water communities 6 . Here, we quantify the depth of the present-day Southern Ocean aragonite saturation horizon using hydrographic and ocean carbon chemistry observations, and use a large ensemble of simulations from the Community Earth System Model (CESM) 8,9 to track its evolution. A new, shallow aragonite saturation horizon emerges in many Southern Ocean locations between now and the end of the century. While all ensemble members capture the emergence, internal climate variability may affect the year of emergence; thus, its detection may have been overlooked by ensemble average analysis in the past. The emergence of the new horizon is driven by the slow accumulation of anthropogenic CO 2 in the Southern Ocean thermocline, where the carbonate ion concentration exhibits a local minimum and approaches undersaturation. The new horizon is also apparent under an emission-stabilizing scenario indicating an inevitable, sudden decrease in the volume of suitable habitat for aragonitic organisms.

Rising atmospheric carbon dioxide (CO 2 ) levels resulting from the burning of fossil fuel and industrial and agricultural activities have been abated by CO 2 uptake by the ocean, which has absorbed nearly one-third of the total anthropogenic carbon added to the atmosphere [10][11][12] . As the ocean absorbs atmospheric CO 2 , its pH and carbonate ion concentration

decrease, thereby decreasing the saturation state Ω = ∕ + -

, which is the product of the concentrations of Ca 2+ and CO 3 2-, divided by the solubility product constant (K sp ) for calcium carbonate (CaCO 3 ) minerals, either aragonite (Ar) or calcite (Ca). Ω Ar and Ω Ca are defined as the ratio of the concentration of dissolved carbonate ions in a given solution to the concentration of dissolved ions in a saturated solution of aragonite and calcite, respectively. Aragonite and calcite are thermodynamically favoured to dissolve once Ω falls below the thermodynamic threshold Ω = 1 and the depth at which this happens within the water column is referred to as the saturation horizon. Ocean acidification makes it harder for marine calcifying organisms (for example pteropods, corals, coccolithophores or foraminifera) to form and maintain their shells 1,7,13 . While pteropods exhibit a physiological negative response between Ω Ar = 0.94 and Ω Ar = 1.12 (ref. 7 ), soft clams, for example, are sensitive to a decrease in Ω Ar well above this thermodynamic threshold 14 .

The Southern Ocean, defined as the region stretching from the Antarctic coastline to 40° S, is especially vulnerable to the effects of acidification relative to lower latitudes. Here, colder temperatures enhance the solubility of CO 2 and persistent upwelling brings carbon-rich water to the surface ocean 1,3,15 . With current CO 2 emission rates, models project that the Southern Ocean's surface will be undersaturated with respect to aragonite by the end of the twenty-first century 1,2,10 . This indicates that key marine calcifying organisms, such as those listed above, may not be able to cope well with future environmental conditions, which could change food web dynamics and have cascading effects on global ocean ecosystems 3,13,15 . Ecosystem impacts in the Southern Ocean will serve as a bellwether for prospective impacts at mid and low latitudes where ocean acidification is projected to occur more slowly 3 .

Here, we use annual output from the CESM Large Ensemble (CESM-LE) 8,9 to study the evolution of the aragonite saturation state under the high-emission Representative Concentration Pathway 8.5 (RCP8.5) (ref. 16 ) scenario (see Methods). The CESM is a stateof-the-art coupled climate model that simulates a unique climate trajectory in each ensemble member 8 . The large ensemble enables a robust estimate of the model's forced response to a given emission scenario and an evaluation of the spread in the response due to internal variability. We focus on the change in the saturation state of the CaCO 3 mineral aragonite, since it is more soluble than calcite at all temperatures and pressures in the ocean and will reach undersaturation earlier.

The depth of the present-day (defined throughout this work as year 2002) observed Southern Ocean aragonite saturation horizon exceeds 1,000 m across most of the basin. Within the core of the Antarctic Circumpolar Current (ACC), we find shallower saturation horizons (~400 m; Fig. 1a). The upwelling of deep water, which contains high CO 2 concentrations from remineralized organic matter, leads to elevated concentrations of dissolved inorganic carbon (DIC) and establishes a naturally shallow saturation horizon in the core of the ACC (refs. 17,18 ). The deepest aragonite saturation horizon depths (~1,400 m) occur in the southwestern Indian Ocean, northeast of coastal Argentina and east of New Zealand.

CESM-LE exhibits a deeper present-day aragonite saturation horizon than that identified by the hydrographic and ocean carbon Letters NATurE CliMATE CHANgE chemistry observations 19,20 (average bias 522 m; Supplementary Fig. 1). To correct for this bias, we use a procedure that pins the model projections to present-day observed distributions of carbonate chemistry, nutrients, temperature and salinity (see Methods). Hereafter, we refer to the bias-corrected model output. This bias correction procedure has been used in the past with much success 1,11 . Moreover, it allows us to describe changes in the saturation horizon due to changes in DIC alone.

The CESM-LE ensemble mean depth of the aragonite saturation horizon, in the locations of the Southern Ocean (south 40° S) where present-day hydrographic data are available, is 83 m in 2100 (Fig. 1c), conforming to results of other recent studies 1,2,10 . Annual average surface ocean aragonite undersaturation begins as early as 2006 in a few discrete locations. Aragonite undersaturation is projected across ~20% of the Southern Ocean surface by 2060, across ~60% of the surface by 2080 and >80% of the surface by 2100.

The CESM-LE projects the emergence of a new shallow saturation horizon across many locations in the Southern Ocean. This emergence is indicated by a step-change in saturation horizon depth of 400 m yr -1 or greater. In some locations, a step-change of as much as 1,000 m in a single year (Fig. 2) is projected. The depth and year of emergence varies spatially, reflecting both natural variation in the present-day saturation horizon depth and spatial variability in the physical circulation of the Southern Ocean. In the core of the ACC in the South Atlantic, we observe the largest step-changes in saturation horizon, ranging from 400 to 1,000 m yr -1 (Supplementary Fig. 2). The step-change is more moderate in the Indian sector, with the exception of a few points near the sea ice edge at 82.5° E.

Step-changes of 500 m yr -1 or more are found throughout the Pacific Sector, extending into the subtropical latitudes. Step-changes also occur for a range of Ω Ar thresholds (Supplementary Fig. 5).

The year of emergence of a shallow aragonite saturation horizon can vary across ensemble members, owing to their different representations of internal variability (Fig. 2, Supplementary Fig. 2), such as the El Niño Southern Oscillation and the Southern Annular Mode which can affect surface

(refs. 21,22 ). For example, Fig. 2a illustrates that while all ensemble members project the emergence of a shallow saturation horizon at 0.5° E and 52.5° S, the year of emergence occurs as early as 2006 in one ensemble member and as late as 2038 in another. This internally-driven spread in the year of emergence means that the average change in the saturation horizon (the mean across all ensemble members) is more moderate at this location. Similar conclusions can be drawn at other locations (Fig. 2b-e), suggesting that using the ensemble mean of several projections from one or more models (as is common practice in the Intergovernmental Panel on Climate Change reports and related publications) may misrepresent the emergence of a shallow horizon and the critical depth where this occurs.

The emergence of a shallow aragonite saturation horizon can be explained by the slow accumulation of anthropogenic carbon in the Southern Ocean thermocline that drives a local reduction of

minimum (Fig. 3). The highest concentrations of

are naturally found in the surface ocean and the lowest concentrations in the bottom of the water column, with a local minimum in the thermocline (Fig. 3c). This (ref. 23 ). In the Southern Ocean, the thermocline minimum in

approaches the saturation concentration for mineral aragonite

, which is primarily a function of pressure and increases with depth in the ocean 23 (Fig. 3c). Thus, an incremental addition of anthropogenic DIC to the thermocline has the potential to lower the

sat(arag) threshold, creating a sudden Ω = 1 horizon in the thermocline. This is illustrated at 0.5° E and 52.5° S, where a small increase in thermocline DIC from 2041 to 2042 causes a new saturation horizon to appear at a depth of 200 m (Fig. 3). Locations that fall within the region impacted by projected sea ice melt (for example, 32.5° E and -65.5° S, Fig. 2d) lack the carbonate ion minimum in the thermocline. Rapid undersaturation of surface waters here is driven by the invasion of anthropogenic DIC and/or by changes in the distribution of natural DIC as rapid ocean warming and freshening affects stratification and ventilation. Because of the technique we used to propagate the bias correction (see Methods),

Observations normalized to year 2002 CESM-LE ensemble mean, 2002 CESM-LE ensemble mean, 2100

internal variability and externally-forced changes in temperature, salinity, alkalinity and nutrients have no direct consequences on the depth of the horizon. However, internally-and externally-driven changes in ocean circulation can affect the interior ocean distribution of DIC and thus indirectly impact the depth of the aragonite saturation horizon.

0 .8 0.9 . 0.9 0.9 0.9 0.9 0.9 0. 9 0.9 0.9 0.9 0.9 0. 9

An ensemble of CESM simulations run under the stabilizingemission scenario RCP4.5 indicates that the emergence of a shallow saturation horizon is unavoidable across a large swath of the Southern Ocean, although the year of emergence can be delayed substantially (Fig. 4). This medium ensemble (CESM-ME, socalled because it has nine ensemble members, see Methods) simulates a similar range of internal variability in the depth of the saturation horizon but with a slower increase in anthropogenic DIC in the Southern Ocean thermocline than that of CESM-LE (RCP8.5). The emergence of a shallow aragonite saturation horizon (defined as the first year where a step-change of saturation horizon is greater than 500 m yr -1 ) occurs approximately 20 years later in CESM-ME (RCP4.5) compared to CESM-LE (RCP8.5). Nevertheless, increases in thermocline DIC occur throughout the southern-most South Atlantic and Indian basins, causing the emergence of a shallow horizon in all CESM-ME (RCP4.5) ensemble members (Supplementary Fig. 3). Across the subtropical South Pacific, where the emergence of shallow saturation horizons were projected in all of the CESM-LE (RCP8.5) ensemble members by 2080 (Supplementary Fig. 2), the CESM-ME (RCP4.5) shows no emergence of a shallow horizon (Supplementary Fig. 3), probably because it occurs later than 2080 (which is the end-date for CESM-ME simulations).

Our analysis implies that Southern Ocean acidification-sensitive organisms will experience a sudden decrease in the volume of their suitable habitat, including shelled pteropods 1,5-7 , foraminifers, cold-water corals 3,24 , sea urchins, molluscs 3 and coralline algae 1,3,24 . Shelled pteropods, the major planktonic producers of aragonite, might be especially vulnerable to these changing conditions since they typically live in the upper 300 m and form an integral component of polar and subpolar food webs [5][6][7] . Pteropods account for a large portion of the flux of CaCO 3 to the deep ocean in the Southern Ocean 25,26 and therefore a decrease in pteropod populations would decrease the amount of CaCO 3 (and, thus, alkalinity) exported to depth. Increased alkalinity remaining in the upper ocean could allow increased oceanic absorption of atmospheric CO 2 , an important negative feedback on climate change. Due to the rapid progression of ocean acidification, pteropods may have a limited time to adapt to a corrosive environment since they produce only two generations per year 27 . While the emergence of a shallow saturation horizon has been projected in coastal upwelling systems 28 , the Southern Ocean is characterized by much lower natural variability in surface ocean 21,29 ). Given this low background variability, organisms in the Southern Ocean may not be able to contend with sudden changes in the volume of their habitat, with far-reaching consequences for fisheries, economies, and livelihoods.

Due to the lack of ship-board wintertime observations, the CESM aragonite saturation horizon is unable to be verified during winter months. Therefore, this analysis focuses only on the annual mean values of aragonite saturation state in the Southern Ocean. Other studies 4,30 , however, show an intense surface wintertime minimum in -CO 3 2 south of the Antarctic Polar Front, which, combined with increasing amounts of anthropogenic CO 2 , will probably lead to earlier undersaturation events during winter. Finally, we note that while CESM-LE and CESM-ME do not represent the potential physiological responses of organisms to ocean acidification, such as altered calcification rates, N 2 fixation, and net primary production, these may also cause future changes in local carbonate chemistry with potentially important climate-carbon feedbacks 15 .

Any methods, additional references, Nature Research reporting summaries, source data, statements of data availability and associated accession codes are available at https://doi.org/10.1038/ s41558-019-0418-8.

Hydrographic and carbon chemistry observations. We use global mapped climatologies of ocean biogeochemical and physical variables collected via hydrographic cruises to identify the present-day Southern Ocean aragonite saturation horizon. DIC and alkalinity are taken from an adaptation of the Global Ocean Data Analysis Product for Carbon, version 2 (GLODAPv2) mapped product 22 that excludes artificial data along the GLODAPv2 mapping boundary at 20° E and includes only data that were quality controlled (that is, no profiles with a maximum sampling depth shallower than 1500 m and no profiles without crossovers) 19 . DIC observations were normalized to the year 2002 before mapping, by removing the temporal trends in DIC and pH due to anthropogenic influence 19,20 . We used mapped climatologies of temperature, salinity, silicate, and phosphate from the WOA 2009 (refs. [31][32][33] ). GLODAP and WOA mapped products are on 1° ×1° grids with 33 standard depth surfaces but here we only used the values in locations where there are observations. We used Mocsy 1.0 (ref. 34 ), a Fortran 90 package that determines the ocean carbonate system, to compute the annual mean saturation state of aragonite at every location and depth in the Southern Ocean. were carried out to 2100, while CESM-ME (RCP4.5) simulations were carried out to 2080. Each ensemble member has a unique climate trajectory because of small round-off level differences in their atmospheric initial conditions 8 . All the CESM ensemble members began with an 1850 control simulation with constant preindustrial forcing. The ocean model physical state was initialized to observations, while the ocean biogeochemical fields were initialized to a state derived from a separate 600-year spin-up. While these spin-ups resulted in a quasi-equilibrium for ocean biogeochemistry, we found significant biases in modelled, present-day Southern Ocean DIC as compared to observations (Supplementary Fig. 1; see also Long et al. 39 ). We therefore used the procedure outlined in Orr et al. 1 and Ciais et al. 11 to make bias-corrected projections of the Southern Ocean aragonite saturation horizon from the two CESM ensembles. For each ensemble member and each projection year, we interpolated the model output to the GLODAP grid and calculated the annual mean DIC anomaly relative to the model estimate in 2002. We propagate this bias correction to 2100 in each ensemble member by adding the simulated model perturbations of DIC, relative to 2002, to the GLODAPv2 DIC climatology, while holding alkalinity, nutrients, temperature and salinity constant. As for the observations (see above), we used Mocsy 34 to calculate the resulting Southern Ocean aragonite saturation state from the bias-corrected DIC model projections. Here too the alkalinity, temperature, salinity, silicate and phosphate were all held constant at their present-day climatological values.

For a given year and desired depth level, ensemble mean values of the simulated variables were computed by averaging across ensemble members. Areas that on an annual average are covered in sea ice were omitted from our analysis due to well-known biases in the present-day CESM sea ice distribution 40 . We define sea ice extent as the northern-most grid point where the simulated sea ice fraction either equals or exceeds 0.2. Seasonal bias. Since the Southern Ocean, due to its remoteness and prohibitive wintertime weather, is almost exclusively sampled during austral summer (December-March), the ship-based biogeochemical observations in GLODAPv2 contain a seasonal bias and very few grid points have data from all seasons. Even when data are available from all seasons, they are often collected many years apart, and these inter-annual variations challenge our ability to identify true seasonal variability. Despite studies showing that seasonal variations of temperature, surface mixed layer depth, and spring blooms have a noticeable impact on Ω Ar and Ω Ca in some regions of the global oceans 41 , no attempt has been made to correct for this seasonal bias in the GLODAPv2 mapped climatologies. This is due both to limited data coverage and that such corrections would have to rely on relationships with ancillary variables and different temporal gap-filling methods 20 . The seasonal measurement bias remains one of the largest sources of unquantified uncertainty for the Ω Ar and Ω Ca estimates in the GLODAPv2 mapped climatologies.

Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.