A New Ocean State After Nuclear War

Abstract

Nuclear war would produce dire global consequences for humans and our environment. We simulated climate impacts of US‐Russia and India‐Pakistan nuclear wars in an Earth System Model, here, we report on the ocean impacts. Like volcanic eruptions and large forest fires, firestorms from nuclear war would transport light‐blocking aerosols to the stratosphere, resulting in global cooling. The ocean responds over two timescales: a rapid cooling event and a long recovery, indicating a hysteresis response of the ocean to global cooling. Surface cooling drives sea ice expansion, enhanced meridional overturning, and intensified ocean vertical mixing that is expanded, deeper, and longer lasting. Phytoplankton production and community structure are highly modified by perturbations to light, temperature, and nutrients, resulting in initial decimation of production, especially at high latitudes. A new physical and biogeochemical ocean state results, characterized by shallower pycnoclines, thermoclines, and nutriclines, ventilated deep water masses, and thicker Arctic sea ice. Persistent changes in nutrient limitation drive a shift in phytoplankton community structure, resulting in increased diatom populations, which in turn increase iron scavenging and iron limitation, especially at high latitudes. In the largest US‐Russia scenario (150 Tg), ocean recovery is likely on the order of decades at the surface and hundreds of years at depth, while changes to Arctic sea‐ice will likely last thousands of years, effectively a “Nuclear Little Ice Age.” Marine ecosystems would be highly disrupted by both the initial perturbation and in the new ocean state, resulting in long‐term, global impacts to ecosystem services such as fisheries.

Introduction

There are more than 13,000 nuclear weapons in the world controlled by nine nations (SPIRI, 2020). States with nuclear weapons, including the United States, Russia, China, India, and Pakistan, have recently embarked on plans to modernize or expand their nuclear arsenals, while North Korea has developed, and Iran is at risk of developing new nuclear weapons capabilities (SPIRI, 2020). The presence of these weapons creates a risk that they will be launched intentionally or by mistake, by unstable leaders, hackers, or computer failure (Ellsberg, 2017; Perry & Collina, 2020). In addition to the devastating loss of human life (Robock et al., 2019; Toon et al., 2019), it is important to quantify the expected collateral damage to the Earth system and all its inhabitants of such a nuclear crisis. Here, we examine the hysteresis response of the ocean to global cooling events driven by the detonation of nuclear warheads.

Hysteresis is the property of a system that has a response that “lags behind” the forcing variable, meaning that the system depends on the history of the forcing (e.g., Visintin, 2013), either resulting in a new state when forcing is removed (irreversible or rate-independent hysteresis), or slowly returns to the initial state (dynamic or rate-dependent hysteresis). In Earth system science, hysteresis has been discussed in reference to the global ocean overturning circulation (Kageyama et al., 2010; Rahmstorf, 1995) and the Antarctic ice sheet (Garbe et al., 2020; Pollard & DeConto, 2005), which both have multiple steady states depending on the forcing history. Other studies have defined rate-dependent hysteresis as not returning to the unperturbed state within a time frame relevant for humans or other organisms (Jeltsch-Thömmes et al., 2020), as in the asymmetry of the climate response to increasing then decreasing CO₂, a temporal lag largely caused by the slow adjustment timescales of the ocean physical and biogeochemical state (Hofmann et al., 2019; Jeltsch-Thömmes et al., 2020; Steffen et al., 2018; Yang & Zhu, 2011). Both irreversible and dynamic hysteresis occur in terrestrial (Drüke et al., 2021) and marine ecosystems under climate disturbances, where they are often termed regime shifts, sometimes due to “tipping points” in the system response to forcing (Brierley & Kingsford, 2009; Cooper et al., 2020; Cowan et al., 2008; deYoung et al., 2008; Drüke et al., 2021). Here, we investigate the hysteresis of the ocean physical, biogeochemical, and lower trophic level ecosystem response to nuclear cooling events.

In the Earth's past, abrupt global cooling driven by volcanic eruptions has triggered atmospheric perturbations that persist for several decades to hundreds of years or longer (Otto-Bliesner et al., 2015; Sigl et al., 2015), with recorded impacts on human civilization including massive famine and fall of empires (Oppenheimer, 2011; White, 2011). Despite the ocean's large role in climate (Rahmstorf, 2002; Sarmiento & Gruber, 2006) and food security (FAO, 2018; Kent, 1997), the global ocean impacts of these cooling events, and in particular the ocean biogeochemical and ecosystem response, remain poorly understood. Like volcanic eruptions (Robock, 2000) and large forest fires (Khaykin et al., 2020; Peterson et al., 2021; Yu et al., 2019, 2021), urban firestorms generated during nuclear war using modern arsenals are expected to loft particles into the upper troposphere and lower stratosphere. However, the smoke would have at least a three times longer residence time than volcanic aerosols, leading to more extended radiation anomalies (Coupe et al., 2019; Mills et al., 2008, 2014; Otto-Bliesner et al., 2015; Robock et al., 2007; Sigl et al., 2015; Toon et al., 2019). This long atmospheric residence time of soot is due to its ability to self-loft to high altitude due to solar heating, as observed in recent forest fires (Khaykin et al., 2020; Peterson et al., 2021; Yu et al., 2019).

In recent work (Toon et al., 2019), we simulated a series of nuclear conflicts between India and Pakistan of varying sizes and impacts on global climate and compared them to a larger war between the US and Russia (Coupe et al., 2019; Robock et al., 2007). As with volcanic-driven events, the reduction of sunlight after nuclear war would lead to global cooling, in turn driving a wide range of global changes, altering the living conditions on Earth. Examples include decreased agricultural and fisheries productivity, with large implications for global food security (Jägermeyr et al., 2020; Scherrer et al., 2020; Toon et al., 2019; Xia et al., 2015); increased surface ocean pH and decreased aragonite saturation (Lovenduski et al., 2020); a 7 year El Niño like event (Coupe et al., 2021); reduced ozone and increased UV radiation at the surface (Bardeen et al., 2021; Mills et al., 2008) and temporary expansion of sea ice (Mills et al., 2014). In simulations, global ocean surface cooling penetrated to depth over decades, indicating the potential for long-term impacts from relatively small global cooling events on the physical ocean state (Mills et al., 2014). Further, in the largest war scenario between the US and Russia, we found that after an initial decline in marine net primary production (NPP), NPP eventually became significantly elevated above pre-war levels long after the solar radi returned to normal. These findings suggest that the extreme cooling caused by a nuclear war might have the ability to trigger changes in the ocean system that persist well beyond the initial perturbation.

While it is known that extreme global cooling events can trigger a rate-dependent hysteresis in the ocean physical state (Bardeen et al., 2017; Mills et al., 2014; Slawinska & Robock, 2018), here, we explore this in detail, and report for the first time, that the ocean also enters a new biogeochemical and ecosystem state. Our study aims to investigate what is driving the new state in global marine productivity, how long this new state is likely to persist, if there is something fundamentally distinct about a large war relative to the smaller cooling events, and if other aspects of the ocean are also affected differently in the more extreme cooling case. To answer these questions, we focus on the ocean impacts of a US-Russia nuclear war that deposits 150 Tg of sunlight absorbing black carbon into the upper atmosphere, triggering nuclear winter, simulated in the Community Earth System Model (CESM) v. 1.3 (Coupe et al., 2019). Smaller India-Pakistan regional war scenarios from Toon et al. (2019) are also studied, where soot injection ranges from 5 to 47 Tg, accounting for varying arsenal sizes and weapons yields.

Results Physical Ocean and Climate Response

Regardless of war location and magnitude, soot reaching the stratosphere from nuclear war firestorms is quickly dispersed globally, absorbing sunlight and reducing the solar radiation to the ocean surface. Soot in the stratosphere dissipates slowly, such that solar radiation reaching the surface is reduced for a decade (Coupe et al., 2019; Toon et al., 2019; Figure 1, Figure S1 in Supporting Information S1). In the US-Russia (150 Tg) scenario, shortwave radiation is reduced by 70% (Figure S1 in Supporting Information S1) and the global average surface temperature (including land) declines by 7°C in the first months, reaching a peak anomaly of −10°C in the third year post-war, with more extreme cooling in the Northern Hemisphere (Coupe et al., 2019). Smaller wars result in proportionally less atmospheric cooling over similar timescales (Coupe et al., 2021; Toon et al., 2019).

Figure 1. Hysteresis plot of annual mean global ocean temperature through time in response to radiative forcing anomalies for all nuclear war scenarios, ranging from 5 to 150 Tg of soot delivered to the stratosphere (colors). Shown are (a) global sea surface temperature (SST). (b) Global potential temperature at 100 m depth (T100), below the mean surface mixed layer (see Figure 2). Arrows represent the direction of time. Dots are plotted for each model year for 15 years (2020–2034). Narrow dashed lines for the 150 Tg case (red) indicate 15 additional years in the simulation, showing the model has not returned to the initial point at the end of the simulation. Solar radiation anomalies are highest in the first year, while SST reaches its maximum anomaly in year 3–4, and T100 reaches a maximum anomaly in year 5–6, depending on the forcing level. See Figure S2 in Supporting Information S1 for anomalies and nitrate hysteresis.

The ocean surface responds more slowly to changes in radiation than the atmosphere and land owing to its higher specific heat capacity. The maximum anomaly in sea surface temperature (ΔSST) occurs 3–4 years post war, later for larger war magnitude, with a maximum ΔSST of ∼6°C in the 150 Tg case in Year 4 (Figure 1a, Figure S1b in Supporting Information S1). The maximum anomaly in subsurface ocean potential temperature (∼4°C at 100 m depth) does not occur until Year 6 for the 150 Tg case (Figure 1b), with later occurrence at deeper depths (Figure 2a). The recovery rate to these temperature perturbations is even slower, such that ocean temperature is still anomalous throughout the water column at the end of the simulations (Figure 1b, Figure S9 in Supporting Information S1). By the end of 2034, 15 years after the war, temperature anomalies are greatest at ∼200 m (Figure S9b in Supporting Information S1), almost −2°C in the US-Russia scenario and over −1°C in the 47 Tg India-Pakistan war. The asymmetry in the initial ocean temperature response time and the long recovery period is a rate-dependent hysteresis, effectively creating a new ocean state for the lifetime of many organisms, including humans, long after the nuclear war-induced radiation anomaly ends. In the discussion we present evidence that this new temperature state is likely to persist for decades near the surface and hundreds of years at depth or longer.

Figure 2. Global ocean state evolution before, during and after a US-Russia (150 Tg) nuclear war. Vertical lines in panels (a and c) indicate the end of the averaging period for profiles (b and d): the war date in 2020 (dotted), mid-cooling event in 2023 (dashed), and 20 years post-war in 2040 (solid). (a) Monthly temperature over the upper water column (colormap) and mixed layer depth (black line); mixed layers increase from 75 to over 200 m in depth during nuclear winter, mixing cold water to the surface and largely eliminating surface stratification. (b) Annual mean temperature profiles; in 2023 (dashed) the upper ocean temperature and thermocline is greatly diminished relative to the pre-war state (dotted); in 2040 stratification has been restored but ocean temperatures have decreased 1°C–2°C throughout the upper water column (solid). Panel (c) same as in panel (a), but for the macronutrient nitrate. (c) Annual mean nitrate profiles; in 2040 (solid) nitrate is elevated throughout the upper water column. Anomaly profiles for temperature and nitrate for all war scenarios are provided in Figure S9 in Supporting Information S1.

Ocean cooling magnitude depends linearly on war severity, both at the surface and below the mixed layer. The magnitude of the initial maximum globally averaged SST anomaly is linearly proportional to the maximum nuclear war driven radiation reduction (Figure 3a), with a −0.55°C anomaly in the annual mean SST per 10 W m⁻² shortwave radiation reduction, where the maximum radiative anomaly is taken from the annual mean timeseries. This climate sensitivity is similar to volcanic forcing response in observations and models (cf., Figure S2 in Chikamoto et al., 2016; Sigl et al., 2015). The range of nuclear war scenarios studied results in −11 to −115 W m⁻² shortwave radiation and −0.5°C to −6.4°C global SST anomalies, with the US-Russia case being the largest. Cooling below the ocean surface mixed layer after radiation anomalies end is also linear with maximum radiation reduction (Figure 3b), with a rate of 0.011°C W⁻¹m² potential temperature perturbation in 2034, 15 years post-war. Regionally, SST anomalies are much more extreme, colder near the coasts and in the Northern Hemisphere due to larger land masses, which cool off faster than the ocean. In the 150 Tg case, maximum SST anomalies exceed −10°C in the West and North Atlantic, and −25°C in the South China Sea (Figure 3c).

Figure 3. Global ocean temperature changes during and after nuclear war driven cooling. (a) Maximum annual global mean sea surface temperature anomaly (ΔSST) versus maximum radiative forcing anomaly for each war simulation. ΔSST decreases linearly (dashed line), at a rate of 0.055°C W−1m2. (b) The global annual mean ocean temperature anomaly at 100 m (ΔT100) in 2034 (Year 15), indicating subsurface cooling in the new ocean state, also decreases linearly with radiation reduction at a rate of 0.011°C W−1m2 (dashed). (c) Maximum monthly ΔSST in the US-Russia (150 Tg) scenario at each model grid cell; ocean temperature anomalies are greatest in the Northern Hemisphere coastal regions and shallow seas. Anomalies are relative to the control mean.

The depth that the surface ocean cooling signal penetrates is modulated by the ocean mixed layer response. Deep vertical mixing at the ocean surface in winter is one of the primary mechanisms facilitating air-sea fluxes and redistributing physical and biogeochemical properties of the upper ocean. Anomalous surface cooling after nuclear war expands the locations and seasons that experience deep mixing (Figure S3 in Supporting Information S1). Normally prevalent for a few months at high latitudes in winter, vertical mixing is deeper and occurs over broader regions during the decadal cooling event, with magnitudes proportional to the forcing (Figure S1b in Supporting Information S1). In the US-Russia case, globally averaged mixed layer depths more than double by 2023 (Figures 2b and 2d); in the Subpolar North Atlantic, winter mixed layer depths triple, and there is no restratification in summer for 5 years (Figures S4b and S4d in Supporting Information S1); similar anomalies occur in the North Pacific and Southern Ocean, extending well into the subtropics in all ocean basins, with less extreme anomalies in the tropics (Figure S3 in Supporting Information S1). The permanent thermocline and pycnocline are eroded due to this enhanced vertical mixing, reducing stratification and thus stability of the water column (Figures 2a and 2b, Figures S4a and S4b in Supporting Information S1).

The increased area and maximum depth of deep mixing drives an exceptional perturbation to the ocean's meridional overturning circulation (MOC). Surface cooling causes a deepened and greatly strengthened Atlantic MOC (AMOC) the first decade after nuclear war (Figure 4, Figure S5 in Supporting Information S1). In the US-Russia case, the AMOC strength at 26.5°N is five times stronger in 2027 than in the control case, well outside the observed range of annual variability (Figure 4a). During the period of greatest response, the level of the maximum AMOC streamfunction is 600 m deeper than in the control: Antarctic Bottom Water (AABW) in the abyssal North Atlantic is displaced by newly formed North Atlantic Bottom Water. The Southern Hemisphere overturning circulation is equally affected by surface cooling: there is increased AABW formation, leading to a temporarily strengthened Southern Hemisphere MOC cell that almost mirrors the Northern Hemisphere's AMOC anomaly across the equator (Figure 4b, Figure S5 in Supporting Information S1). Depending on the latitude, MOC strength peaks 5–10 years after the nuclear war before returning to its pre-war strength. In both hemispheres this intense MOC pulse scales linearly with the amount of soot injected into the atmosphere in each of the scenarios (Figure S6 in Supporting Information S1). In both hemispheres, the MOC streamfunction anomaly is almost entirely due to the resolved ocean transport with only a very minor contribution from the mesoscale and submesoscale parameterizations in CESM (not shown, Fox-Kemper et al., 2011; Gent & McWilliams, 1990).

Figure 4. Evolution of oceanic meridional overturning circulation (MOC) post-war. (a) Atlantic MOC strength at 26.5°N post-war (solid black line) as compared to control simulations (black dashed lines) and observed interannual range of overturning streamfunction from the RAPID-MOCHA time series (gray shading; Smeed et al., 2019). (b) 2025–2030 global (all ocean basins) MOC anomaly (shading, positive clockwise) from the control mean (contours, solid indicates clockwise).

Arctic sea ice grows rapidly in response to the reduced solar forcing (Figure 5). The normal range of sea ice extent and thickness is taken from the CESM1 Large Ensemble (CESM1-LE) preindustrial control, a set of simulations that allow the quantification of the Earth system's range of internal variability (Deser et al., 2020; Kay et al., 2014). In the Northern Hemisphere, sea ice volume grows beyond the range of this internal variability in all but the smallest war scenario, in most cases exceeding the internal variability during the first year post-war (Figures 5d and 5h). In the US-Russia scenario, peak Arctic sea ice extent expands by 10 million km², covering over 50% more area, including normally ice-free coastal regions important for fishing, aquaculture, and shipping across the Northern Hemisphere (Figures 5a–5d). In this scenario, ice expands most dramatically in the Sea of Japan, South China, Caspian, and Baltic Seas, where concentrations of sea ice above 50% are present in normally ice-free regions; sea ice thickness more than doubles Arctic-wide, in many basins deepening more than 2 m (Figures 5e–5h). After the extreme forcing event ends, Arctic sea ice extent and thickness settle into a new steady state in the US-Russia scenario so that the expanded ice is retained until the end of the simulation, outside the range of internal variability (Figures 5d and 5h), indicating a rate-independent hysteresis response. It is likely that the new Arctic state is amplified by the surface albedo feedback, which would lower the summer ice melt (Miller et al., 2012; Perovich & Polashenski, 2012). Southern Hemisphere sea ice remains within the range of internal variability, such that maximum sea ice extent and volume are never larger than in the CESM-LE (Figure S7 in Supporting Information S1).

Figure 5. Post-war Arctic sea ice evolution. Arctic 2020–2025 mean sea ice concentration (%) for (a) the US-Russia Nuclear War (NW) scenario, (b) the control scenario, and (c) the difference in concentration between the two scenarios, and Arctic mean sea ice thickness (m) for (e) the US/Russia Nuclear War scenario, (f) the control scenario, and (g) the difference in thickness between the two scenarios. The Northern Hemisphere annual mean time series of (d) sea ice extent and (h) sea ice volume is shown for all war scenarios (colors) and control scenarios (black), where the Community Earth System Model-Large Ensemble experiment mean (solid gray line) and standard deviation (dashed) over the preindustrial period are given to demonstrate the natural, internal variability within the model (see Section 5).

In summary, although solar radiation and mixed layer depths have returned to their pre-war levels within ∼10 years (Figures S1 and S2 in Supporting Information S1), the ocean physical state remains significantly altered at the end of the simulations (15–30 years), with a cooler water column, altered deep stratification, and increased Arctic sea ice (Figures 1b, 2a, 2b and 5, Figure S4a and S4b in Supporting Information S1). The changes are the most pronounced in the US-Russia case. For example, the 12°C global isotherm is ventilated during the extreme cooling event, shoaling from 200 m to the surface, and has deepened to only 160 m by the end of the simulation, 30 years after the war, leaving a ∼1°C cooling signal throughout the upper water column (Figures 2a and 2b). In the new state, temperature anomalies below the thermocline are linear with the radiation reduction, such that the 100 m global temperature average 14 years after the war is 0.1°C cooler per 10 W/m² (Figure 3b). Changes in the ocean biogeochemical and ecosystem state stem from these physical changes in the ocean.

Biogeochemical and Phytoplankton Impacts

Like temperature, the vertical distribution of biogeochemical tracers is strongly modified by the enhanced vertical mixing and overturning initiated by the extreme forcing event (Figure 2, Figures S8–S9 in Supporting Information S1). In the ocean, nitrate is a primary nutrient controlling ecosystem productivity. In the US-Russia scenario, the global 12 mmol/m³ nitrate isopleth shoals from 125 to 80 m during the first 4 years following the war, then stabilizes at ∼100 m until the end of the simulation, with a very slow recovery rate (<1 m/yr; Figure 2d). This shoaling brings high nutrient water closer to the surface where it can be utilized by phytoplankton. After nuclear cooling ends, nitrate anomalies are greatest below the seasonal mixed layer (Figure 2) with the magnitude proportional to the radiative forcing (Figure S9 in Supporting Information S1). The nitracline in the Subpolar North Atlantic is largely destroyed during the cooling event, and the AMOC intensification drives the transport of low nitrate surface water into North Atlantic Bottom Water (Figure S9 in Supporting Information S1). Low oxygen zones in the ocean have negative impacts on many marine species, and subsurface oxygen is especially limited in the Pacific where water is poorly ventilated (e.g., Kwiatkowski et al., 2020; Schlitzer, 2000). After nuclear winter, high latitude waters in both the North and South Pacific are ventilated past 1,000 m, eliminating a large expanse of the oxygen minimum zone in the North Pacific (Figure S8 in Supporting Information S1).

Photosynthetically available radiation (PAR) is solar radiation to the ocean's surface used for marine photosynthesis (Figure 6a), in CESM equivalent to 45% of the shortwave radiation at the top of the water column (Figure S1a in Supporting Information S1). The global NPP in the ocean (i.e., the new growth of marine algae, which makes up the base of the marine food web) declines quickly after the war in all scenarios, with the magnitude and rate of NPP decline proportional to the reduction in PAR (Figure 6; Figure S1d in Supporting Information S1). In the US-Russia scenario, global NPP is reduced by nearly 50% in the months after the war and remains depressed by 20%–40% for over 4 years (Figure 6a), with the largest reductions in the North Atlantic and North Pacific (Figure S10 in Supporting Information S1). This initial productivity reduction is accompanied by a global shift in phytoplankton community abundance to favor diatoms (Figure S11 in Supporting Information S1), which generally outcompete other plankton types in low-light, high-nutrient conditions (Miller & Wheeler, 2012). By 2040, global NPP has recovered and stabilized at an elevated state in the US-Russia scenario, ∼6% higher than before the war (Figure 1), with gains as high as 20%–30% across the tropics and subtropics (Figure S10 in Supporting Information S1). The regional NPP response to anomalous forcing is consistent across all scenarios, with the magnitude of the elevated NPP in this new state again proportional to the degree of anomalous radiation forcing (see Scherrer et al., 2020). The magnitude and larger spatial extent of the positive NPP anomalies in low latitudes drive the elevation of globally averaged NPP in the US-Russia scenario, while smaller war scenarios have compensating negative NPP anomalies at higher latitudes and in the Equatorial Pacific.

Figure 6. Change in marine net primary production (NPP) and its drivers, photosynthetically available radiation (PAR) and nutrient availability (nitrate), after a US-Russia (150 Tg) nuclear war. Thin lines in panels a–c are monthly values, thick are annual means, control simulation values are dashed, and a vertical dashed line indicates the war timing. (a) Global NPP decline (black) initially follows PAR reduction (red); after the smoke clears, NPP is elevated by ∼6% due an increase in surface nitrate availability (blue). (b) NPP in the Subpolar North Atlantic (red outline in panels (d and e)) is extinguished (>95% reduction) in the year after the war, takes a number of years to recover, and is not elevated post-recovery. (c) NPP in the Sargasso Sea (blue outline in panels (d and e)) decreases by only 30% in the first year, then is elevated by 20%–30% until the end of the simulation. (d) In years 1–2 post war, NPP reduction is greatest in the Northern Hemisphere and outside of the subtropical gyres. (e) In years 19–20, after the atmospheric forcing has subsided, NPP is elevated across the tropics and subtropics.

Understanding what drives these changes in NPP requires disentangling the various controlling factors of phytoplankton growth. Simulated phytoplankton cell division rates are regulated by light, temperature, and nutrient limitation terms (see Methods). Each limitation term scales the maximum division rate for each of the three simulated phytoplankton functional types, such that lower values translate to slower growth rates. Simulated nutrient limitation is determined by the most limiting nutrient: nitrate, phosphate, or iron for all phytoplankton, and silica for diatoms (Moore et al., 2013). Nuclear war increases light limitation by both reducing PAR and deepening the mixed layer, with the latter transporting phytoplankton out of the photic zone. Cooling increases temperature limitation, both directly and through driving deeper mixing and shoaling of the thermocline. In contrast, nutrient limitation is generally relieved by enhanced vertical mixing and shoaling of the nutricline. Perturbations to the three limitation factors have different timescales and different levels of importance in different ocean biomes and through time.

Globally, light is the dominant limiting factor during the cooling event, while nutrient limitation controls production changes in the new ocean state. During the cooling event, reduced PAR and deeper mixing increase phytoplankton light limitation, which along with increased temperature limitation results in decreased global phytoplankton division rates, biomass and NPP, despite an increase in the vertical nutrient flux and thus relieved nutrient limitation (Figures 2, 6 and 7). In the new ocean state, light limitation is restored to pre-war levels, however the shoaling of the thermocline and nitracline result in elevated temperature limitation and relieved nutrient limitation (Figure 7). In the new ocean state, nutrient relief dominates, resulting in an ocean where globally averaged phytoplankton division rates, biomass stocks, and NPP are all elevated (Figures 6 and 7).

Figure 7. Globally averaged drivers, mechanisms, and response of phytoplankton anomalies. Monthly (thin) and annual means (thick) for the nuclear war case (solid) and control climate (dashed). A vertical dashed line denotes the war date. Time series are plotted for the physical drivers (a) mixed layer depth (MLD) and (b) surface photosynthetically available radiation (PAR). In the middle, time series are plotted for the (c) biomass-weighted limitation terms that control phytoplankton growth: nutrient, temperature, and light limitation; limitation terms are multiplicative and lower numbers are more limiting. Below, time series are plotted for the biological response of (d) phytoplankton division rates and (e) phytoplankton biomass. Initially, light limitation is the primary driver of production decline (Figure 6a), with contribution from temperature limitation. Post-war, reduced nutrient limitation results in higher primary productivity (Figure 6), despite increased temperature limitation.

NPP anomalies during and after the cooling event exhibit strong regionality, caused by differing balances in the relative limitation factors (Figure 6, Figures S11 and S12 in Supporting Information S1). Short-term post-war NPP reductions driven by light limitation are most pronounced at high latitudes, where phytoplankton are already light limited and often constrained to short, seasonal blooms (Racault et al., 2012; Sundby et al., 2016). For the Subpolar North Atlantic in the US-Russia case (Figure 6b, Figures S4 and S12a–S12e in Supporting Information S1), surface PAR is almost totally extinguished immediately post-war and mixed layer depths reach over 1,000 m, leading to intense light limitation and causing NPP to shut down completely for an entire year. The spring bloom is subsequently stifled for another 2 years (by ∼75% and ∼33%, respectively) and occurs later for 5 years (Figure 6b).

In contrast to the Subpolar North Atlantic, NPP in the Sargasso Sea in the US-Russia scenario (Figure 6c, Figures S12f–S12j in Supporting Information S1) decreases much less in the first year of the cooling event (∼30% reduction, ∼10 mmol C m⁻² d⁻¹), with some areas showing increased NPP during the cooling event, after which there are 5 years of summer blooms up to double the size of control levels in 2023. This region eventually settles to a new state that is ∼20% more productive than the control run. Similarly, significant long-term NPP enhancement is evident throughout the subtropical gyres (Figure 6e, Figure S10 in Supporting Information S1). Higher NPP in the new state is driven by partial relief of the persistent nutrient limitation (Figure 7c, Figure S12h in Supporting Information S1; Chavez et al., 2011; Longhurst, 2006). In these biomes limiting nutrients (nitrate and phosphate; Harrison et al., 2018) are delivered by vertical mixing (Menzel & Ryther, 1959) and this nutrient delivery is enhanced by the shoaled nutricline in the new ocean state, stimulating faster phytoplankton division rates, larger biomass stocks, and elevated NPP.

War-driven changes in iron limitation affect productivity differently in different ocean biomes. After light is restored, NPP is generally elevated in the low latitudes, but not at high latitudes or in the Eastern tropical Pacific, despite shoaling of the nitracline there (Figure 6, Figures S4, S9 and S10 in Supporting Information S1), indicating other nutrients play a role. Ocean biomes can be roughly divided into iron limited, such as the high latitude and Eastern Pacific oceans, and nitrate or phosphate limited (Harrison et al., 2018; Misumi et al., 2014; Moore et al., 2013). While nitrate and phosphate biogeochemical cycling and vertical distributions are similar, iron cycling is distinct, complicated by the affinity of iron to bond with marine detritus, which sinks and removes bioavailable iron from the water column (Misumi et al., 2011; Moore et al., 2004).

At high latitudes, and in the Eastern tropical Pacific, we find that an increase in diatom biomass, driven by the shoaling of the nitracline, also increases detritus production and in turn iron scavenging, depleting iron and driving enhanced iron limitation. After the war, both iron and nitrate limitation are relieved globally through a combination of enhanced vertical nutrient delivery and the light-limited decline in production, leaving available nutrients unutilized (Figure 7, Figures S11 and S12 in Supporting Information S1). Nitrate limitation continues to be relieved after the cooling event, following from the shoaling of the nitracline (Figure 2, Figure S9 in Supporting Information S1). In contrast, iron limitation is only reduced for a few years, both regionally and globally, then iron limitation increases to a more limiting state than before the war, remaining more limiting until the end of the simulation (Figures S11b and S12c in Supporting Information S1). The low light and high nutrient levels during the cooling event favor the production of diatoms over small phytoplankton, such that the global diatom biomass is elevated by 50% in 2023 (Figure S11a in Supporting Information S1). This shift to diatoms leads to an increase in particulate organic carbon (POC) production and in turn an increase in iron scavenging, depleting the water column of iron (Figure 8; Moore et al., 2001; Misumi et al., 2011, 2014). In the iron-limited ocean biomes such as the North Atlantic, this increased iron limitation reduces the ability of phytoplankton to fully utilize the enhanced nitrate. In the new ocean state spring, when iron levels are highest, higher nitrate results in larger diatom blooms (Figure S11 in Supporting Information S1). Detritus production is in turn enhanced nonlinearly, as POC production in CESM-biogeochemical elemental cycling (BEC) grows quadratically with plankton biomass (see Harrison et al., 2018). POC bonds with iron, increasing the scavenging rate (Figure 8), and maintaining higher levels of iron limitation in the new ocean state (Figures S11b and S12c in Supporting Information S1). Thus, the shoaled nitracline only translates to increased production in regions which are iron-replete, like the subtropical gyres (Figure 6, Figures S10 and S12 in Supporting Information S1). This new, more iron-limited global ocean state continues until the end of the simulation.

Figure 8. Dissolved iron (Fe) and dissolved iron scavenging rate anomaly vertical profiles in the US-Russia case for the Subpolar North Atlantic (Figure 6) before, during and after nuclear cooling. (a) Iron is higher throughout the upper water column during the cooling event in 2023 (dashed), than pre-war (2019; dotted), due to the enhanced vertical mixing and lack of uptake by phytoplankton production. By 2034 (solid) Fe is significantly lower than the pre-war levels. Iron scavenging rate anomaly (b) through the entire water column and (c) in the upper 500 m. In the year of the war (2020, dotted), scavenging declines due to lack of production, then recovers (2023, dash dot), so that by 2027 (dashed) scavenging is elevated throughout the water column. Twenty years after the war (2040, solid) scavenging remains elevated, especially in the upper water column.

Discussion

How long would it take the ocean physical and biogeochemical state to return to what it was pre-war? We expect the new ocean state not to be a new steady state, that is, an irreversible hysteresis, but an extremely long transient, with a few exceptions. Vertical distributions of physical and biogeochemical tracers take hundreds to thousands of years to achieve steady-state for a given climate state (DeVries & Primeau, 2011; Sarmiento & Gruber, 2006). Ocean temperature and density perturbation adjustment timescales vary by depth and ocean basin, and are a function of the ocean's internal variability and key circulation components such as MOC, which are in turn affected by climate perturbations (Hogan & Sriver, 2019). The ocean response to such perturbations is asymmetric, with a generally faster response to cooling than to warming due to cooling-enhanced convection and overturning (Yang & Zhu, 2011), as seen in our simulations. The precise response to a given perturbation is difficult to predict a priori; there is no fundamental timescale that separates rapid and slow adjustments (e.g., Sherwood et al., 2015), or separates linear transient responses from those that arise due to climate feedbacks (Knutti & Rugenstein, 2015). Estimates for a fully equilibrated oceanic response to an increase in atmospheric CO₂ range from 200 years for the upper ocean (Yang & Zhu, 2011) to thousands of years for the MOC (Jansen et al., 2018), where the latter is essentially a diffusive adjustment timescale of the abyssal ocean. Jansen et al. (2018) note that the adjustment of the MOC in their model occurs on three distinct timescales, where the strength of the MOC is partially restored after only a few decades. Other adjustments, such as to the thermocline density, Southern Ocean eddy-driven circulation (e.g., Allison et al., 2011; Jones et al., 2011), and abyssal circulation rely on longer climate feedbacks that are beyond the scope of this paper. Based on the simulations here, and on these other studies, the complete recovery from a nuclear war will be spatially inhomogeneous and in some areas on the order of thousands of years, if ever.

With the above caveats in mind, we estimated a recovery timescale as a function of depth for the 150 Tg case (Figure 9). At each depth, this was calculated using the global mean temperature anomaly divided by its recovery rate, where both are taken after the nuclear war-driven radiation anomaly has ended and the rate of temperature change has slowed (Figure S13 in Supporting Information S1). The temperature anomaly was averaged over the last 5 years of the long control run simulation (2036–2040) to account for interannual variability in the control run. The recovery rate (Figure 9a) was computed as the time rate of temperature change averaged over the last 10 years of the 150 Tg simulation (2040–2049), less the time rate of temperature change averaged across the last 10 years of the three control runs to correct for model drift. The estimated recovery timescale (Figure 9b) ranges from 15 to 20 years in the top 100 m, increasing to 80 years between 200 and 400 m, and then approaches a singularity as the recovery rate approaches zero at ∼600 m. At this depth, the ocean in the nuclear war simulation is still cooling, relative to a warming drift in the control simulations (Figure 9a), indicating the nuclear war cooling signal is still penetrating the deep ocean 30 years after the war. This decadal scale cooling at depth underscores the long timescales associated with the ocean response to this extreme cooling event.

Figure 9. Estimated (a) global average recovery rate of the temperature anomaly and (b) timescale of recovery as a function of depth for the 150 Tg scenario. Note difference in vertical scale. (a) The recovery rate of the temperature anomaly (dΔT/dt, solid) was estimated by subtracting the drift rate (see Figure S13 in Supporting Information S1) averaged over 10 years of the control mean (dotted) from the rate of temperature change over the last 10 years of the 150 Tg simulation (dashed). (b) The recovery timescale (tR) was estimated using the global mean temperature anomaly divided by the recovery rate relative to the control drift, where both are taken after the nuclear war-driven radiation anomaly has ended. The recovery timescale approaches a singularity at 600 m, as the recovery rate approaches zero.

For Arctic sea ice, the hysteresis does not just result in a long transient but is likely a new steady state for many nuclear cooling scenarios, as in historical global cooling events. While only the US-Russia scenario was simulated long enough to achieve a steady sea-ice state, it is likely that sea ice volume would remain elevated in the four nuclear war simulations above 5 Tg, given their elevated state at the end of the cooling event (Figure 5). Similar elevated sea ice events are observed in the paleoclimate record after volcanic cooling events. For example, strong evidence supports that an expanded Arctic sea ice state during the Little Ice Age (LIA) was triggered by volcanic forcing within the range of the scenarios here. Decadal scale cooling, driven by tropical volcanic eruptions, has been simulated to induce an albedo feedback (Perovich & Polashenski, 2012) that causes centennial scale increases in Arctic sea ice during the LIA (Zhong et al., 2011), a finding supported by a range of proxy records across the Arctic (Miller et al., 2012). Climate proxy records (Sigl et al., 2015) indicate that the reduction in solar forcing of the largest pre-LIA volcanic eruption was 35 W/m², a small fraction of the 120 W/m² reduction in our US-Russia nuclear winter simulation, and similar to the 16 Tg India-Pakistan case (Coupe et al., 2021). Thus, both the short-term and persistent expanded sea ice seen in simulations of nuclear war driven cooling is consistent with historical evidence. We conclude that a large nuclear war would likely induce a nuclear little ice age and be more likely to induce such an event than a volcanic eruption with the same radiation anomaly, due to the shorter lifetime of volcanic aerosols and thus briefer cooling duration (Robock, 2000). Anthropogenic climate change, which was not simulated here, would mitigate some cooling and the resulting ice expansion (Bethke et al., 2017). However, for a large nuclear war anthropogenic CO₂ emissions would likely plummet, eventually slowing or even reversing global warming (Matthews & Weaver, 2010; Robock, 2015).

How would the immediate and long-term impacts of nuclear cooling affect marine ecosystems and fisheries? The large and rapid changes in marine primary production, temperature and other physical conditions would result in declines in marine consumers like fish, shellfish and larger bodied animals, both warm and cold blooded (Guiet et al., 2016; Heneghan et al., 2021; Jennings et al., 2008; Schwartz, 1978). The overall abundance of consumers (i.e., fish) and the presence of large-bodied species is directly related to the magnitude of primary production (e.g., Guiet et al., 2016; Jennings et al., 2008), especially for non-mobile species (Watson et al., 2015), and lower ocean temperatures slow the growth rates of fish (Brown et al., 2004). Our recent modeling using a global fisheries model forced by these simulations found that global fish biomass would decline by ∼20% over the first 10 years after the US-Russia war scenario, proportional to declines in primary production, but also depending strongly on response of fishing pressure after the war (Scherrer et al., 2020). These estimates have high uncertainty (Scherrer et al., 2020), as the fisheries model assumes immediate adaptation of marine communities to rapid climate velocities, and does not resolve species level processes such as how changes in plankton bloom timing would affect feeding and reproduction success of zooplankton and fish, or how light reductions affect visual hunting (Aksnes et al., 2004; Cahill et al., 2013; Durant et al., 2007; Pinsky et al., 2020).

Given the large uncertainties and limited modeling performed for marine consumers, it is highly uncertain what the upper trophic level communities would look like in the new ocean state. Global marine ecosystem models generally simulate functional groups or fish size spectra and not individual species (Heneghan et al., 2021), and the explicit links between the changes in the composition of phytoplankton, zooplankton and fish communities need to be investigated further. While there is evidence that some pelagic fish species rapidly shifted geographic range and adapted to the climatic perturbation following the Tambora volcanic eruption (Alexander et al., 2017), we underline that most of the war scenarios modeled here cause much larger and more sustained climate anomalies. Temperature impacts from nuclear war are particularly extreme in coastal regions where the majority of fish catch and marine ecosystem services are provided (Figure 3c). Benthic and sessile marine organisms in coastal zones would be particularly vulnerable to the locally extreme temperature and NPP shifts (Figures 2 and 6). The long-term thermal tolerance of many species and life history stages would likely be exceeded by the velocity, magnitude and duration of the nuclear war driven climate event, leading to decreased fitness and increased mortality (Peck et al., 2009; Rezende et al., 2011). Thus, the strong declines in temperature and primary productivity projected in many of the nuclear war scenarios could lead to collapses of marine food-chains and potentially extinctions of some endemic marine consumers. The timescale of recovery to such disturbances is unclear.

The immediate and persistent increase in tropical and subtropical primary production after nuclear cooling would have multiple impacts. To offset the inevitable agricultural losses and famine if fisheries collapse following a nuclear conflict, society will need to look to alternative food sources, such as wild caught and farmed seaweed. In this scenario, the immediate and persistent increase in tropical and subtropical primary production after nuclear cooling will likely provide the best region for seaweed production. Building on its usage as a famine food throughout history, using wild caught and farmed seaweed after extreme climate events (Baum et al., 2015; Mouritsen et al., 2021) has been proposed as an alternative food source to offset inevitable agricultural losses and famine following a nuclear conflict (Jägermeyr et al., 2020; Xia et al., 2015). The best region for seaweed production would likely be in areas where primary production is not impacted or is higher due to increased nutrient supply, that is, in many regions across the tropics (Figure 6), and findings here support seaweed production could be accomplished both immediately during nuclear cooling (Figure 6d) and in the recovery phase (Figure 6e). Additionally, a long-term increase in wild seaweed productivity in the subtropics could increase habitat for a great diversity of marine animals (Laffoley et al., 2011; Witherington, 2002), but also act as marine debris, negatively impacting fishing and ship traffic (Partlow & Martinez, 2015; Resiere et al., 2018). Overall, enhanced NPP following the initial cooling is expected to increase the food available for consumers, and modeling suggests about a 15% increase in global fish biomass by year 15 post-war relative to control simulations in the 150 Tg case (Scherrer et al., 2020), keeping in mind the caveats outlined above. However, large blooms in the cooling recovery period could likely include harmful algal blooms, which are toxic to marine life, including birds, fish and corals (Abram et al., 2003; Anderson et al., 2002), as well as humans (Kirkpatrick et al., 2004). The combination of increased productivity in the tropics and increased iron limitation at high latitudes would impact large scale biogeochemical cycling (Krumhardt et al., 2020; Moore et al., 2018; Sarmiento & Gruber, 2006) and there is also indication that less iron availability would impact fishes directly (Galbraith et al., 2019). We can surmise that the number and character of ecosystem disruptions in the new ocean state, driven either through nuclear war or volcanic-driven cooling, would be extensive, persistent, diverse, location specific, and difficult to predict due to nonlinear ecosystem interactions.

Conclusions

Following each of the nuclear war scenarios, a decadal solar radiation reduction and cooling event ensues. The low light and rapid, sustained cooling cause large physical perturbations to the ocean, including an intensification of vertical mixing, enhanced overturning and expansion of sea ice. In all nuclear war scenarios, these perturbations drive a hysteresis, generating a new ocean state, where density stratification and biogeochemical tracers are altered throughout the water column. The magnitude of the new state and timescale of recovery are affected by the magnitude of the cooling anomaly and depend on the system effected. In all nuclear war scenarios, modifications to temperature and biogeochemical profiles persist for many decades, and likely hundreds of years, owing to the long recovery timescales of the deep ocean. Increased Arctic sea ice extent and volume are likely permanent in US-Russia (150 Tg) war simulation, as in the case of the volcanic cooling driven initiation of the Little Ice Age at the end of the thirteenth Century CE. Increased macronutrient availability, driven by changes in biogeochemical profiles, boosts marine productivity in the tropics and subtropics during nuclear cooling, and even more once nuclear cooling ends. However, increased iron limitation at high latitudes due to an altered phytoplankton community results in no enhancement of production there, despite the higher nitrate. Enhanced tropical productivity in the US-Russia scenario leads to an increase in global NPP that persists for decades. These results indicate major changes in ocean ecosystems after a global cooling event, both from the nuclear cooling event itself, but also from long-term changes to the marine physical and biogeochemical system. These findings expand our understanding of how both historical and potential future global cooling events impact ocean physics, biogeochemistry, and ecosystems.

Methods

Simulations used the CESM version 1.3, a state-of-the-art coupled climate model consisting of atmosphere, ocean, land, and sea ice components. Ice is simulated using the Los Alamos sea ice model, CICE, version 4 (Hunke & Lipscomb, 2008) which is used in CESM1 as described by Hurrell et al. (2013). CICE includes simulation of deposition and cycling of aerosols (dust and black carbon) within the ice pack and their resultant impacts on albedo and melting (Holland et al., 2012). CESM1 implements the Parallel Ocean Program (POP) physical ocean model (Danabasoglu, Bates, et al., 2012), here at nominal 1° horizontal resolution and with 60 vertical levels, and the BEC ocean ecosystem-biogeochemistry module, which represents the lower trophic levels of the marine ecosystem, and a dynamic iron cycle (Moore et al., 2004, 2013, 2018; Harrison et al., 2018; Lindsay et al., 2014; Long et al., 2013; Rohr et al., 2017). Similar to other Coupled Model Intercomparison Project v.5 (CMIP5) class models (Laufkötter et al., 2015, 2016), BEC simulates three phytoplankton functional types: diatoms, small phytoplankton, and diazotrophs; the productivity (carbon fixation) of the three phytoplankton groups are combined to generate NPP (Krumhardt et al., 2017). The CESM-BEC ecosystem and biogeochemistry model is well-validated in a variety of scenarios and performs favorably when compared with other CMIP class models (e.g., Harrison et al., 2018; Rohr et al., 2020; Tagliabue et al., 2016 and references therein). AMOC strength has been well-validated in CESM1 simulations (Danabasoglu, Yeager, et al., 2012; Maroon et al., 2018; Yeager & Danabasoglu, 2012) and AMOC strength at 26.5N in control simulations presented here is comparable to the observed AMOC RAPID-MOCHA array (Figure 3; Smeed et al., 2019).

Atmospheric circulation and chemistry are simulated in this version of CESM1 using the Whole Atmosphere Community Climate Model (WACCM; Marsh et al., 2013) with nominal 2° resolution and 66 vertical levels, a model top at ∼145 km, and uses the Rapid Radiative Transfer Model for GCMs (Iacono et al., 2000). The Community Aerosol and Radiation Model for Atmospheres (CARMA; Toon et al., 1988; Bardeen et al., 2008) is coupled with WACCM to simulate the injection, lofting, advection, and removal of soot aerosols in the troposphere and stratosphere, and their subsequent impact on climate (Bardeen et al., 2017, 2021). CARMA has 21 different size bins for the aerosols with different optical properties, allowing particle growth to change the amount of absorption and scattering. Hygroscopic growth is not included, but particles coagulate together at a rate which is partially a function of relative humidity. This provides a realistic simulation of soot growth, transport, and deposition.

The climate response to nuclear war is simulated by injecting black carbon (soot) into the stratosphere in a layer between 100 and 300 hPa over a 1 week period starting on 15 May above the U.S. and Russia, or the South Asian subcontinent (Bardeen et al., 2008; Coupe et al., 2019; Toon et al., 2019) in May in a modern simulation, which here we assume represents the year 2020. In total, six nuclear war scenarios are simulated. Five India-Pakistan nuclear war scenarios are considered which involve soot injections of 5, 16, 27.3, 37, and 46.8 Tg, representing ranges in arsenal sizes. Scenario development, including uncertainties in smoke production and model implementation, are detailed in Toon et al. (2019), along with an overview of global simulation results and direct casualty estimates. These simulations are each run for 15 years. One United States-Russia scenario is considered with a 150 Tg soot injection, for which the simulation was extended to 30 years to assess the ocean recovery to this larger forcing over a longer timescale. This U.S.-Russia scenario is a legacy scenario that assumes both countries use most of their nuclear arsenals on each other's cities (Robock et al., 2007), which is still possible with modern nuclear arsenals. We revisited this scenario in the more modern CESM model and found that the results were very similar to previous simulations (Coupe et al., 2019). An ensemble of three control simulations were run for 19 years to simulate the impacts of interannual variability on physical and biogeochemical processes during the nuclear war driven cooling phase, and to assess the nuclear war climate perturbation relative to any model drift. One of these was extended to 25 years to compare to the longer nuclear war simulation. The nuclear war simulations were branched off the long control run after 4 years.

The ocean was initialized from rest using a standard set of CESM-POP-BEC initial physical and biogeochemical conditions taken from observations and long model spin-ups where observations are not available. The nuclear war scenarios were initiated after 4 years of ocean model spin up, long enough for the ocean physics and productivity response to equilibrate. An ensemble of three control simulations were run to simulate the impacts of interannual variability on physical and biogeochemical processes. Because the ocean simulations were initialized from previously spun-up and observational initial conditions, as is standard for computationally expensive simulations that preclude long-spin up simulations (e.g., Harrison et al., 2018; Liu et al., 2021), model drift was small relative to the nuclear war driven perturbations. Similarly, it was determined that internal climate variability was much less than the nuclear war perturbation, even in the 5 Tg scenario.

In CESM-BEC, phytoplankton net population growth is governed by a photosynthetic NPP term opposed by a loss term. NPP and loss terms are computed independently for each phytoplankton class (diatoms, small phytoplankton, and diazotrophs) and composed of a volumetric specific rate multiplied by the biomass concentration. The volumetric specific division rate is governed by temperature dependence, light availability, and multi-nutrient (N, P, Si, and Fe) limitation (Geider et al., 1998). Multi-nutrient limitation is treated following Liebig's Law of the Minimum (van der Ploeg et al., 1999) such that the maximum specific division rate is only scaled by the most limiting nutrient limitation term. However, light, temperature and nutrient limitations are co-limiting and each multiplicatively impact the maximum specific division rate (see Rohr et al., 2020). All limitation terms vary from 0 to 1 and scale the maximum division rate such that lower values translate to more limitation and reduced growth. The volumetric specific loss rate is composed of a nonlinear grazing rate, a linear mortality rate and a quadratic mortality/aggregation rate. Ocean acidification does not affect phytoplankton NPP or calcification in this configuration of CESM (Lovenduski et al., 2020). NPP and biomass are integrated across all three phytoplankton groups to a depth of 150 m. Division rates and limitation terms are biomass-weighted and depth averaged. Spatially averaged terms are further biomass weighted latitudinally and longitudinally.

Acknowledgments

This work was supported by the Open Philanthropy Project. Data from the RAPID AMOC monitoring project is funded by the Natural Environment Research Council and are freely available from www.rapid.ac.uk/rapidmoc. The authors acknowledge the University of Colorado Boulder Research Computing Group, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University, the Texas Advanced Computing Center (http://www.tacc.utexas.edu) at The University of Texas at Austin for providing visualization resources that have contributed to the research results reported within this paper, as well as high-performance computing support from Cheyenne (https://doi.org/10.5065/D6RX99HX) provided by the NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation. A. DuVivier acknowledges the CESM project, which is supported primarily by the National Science Foundation (NSF) at the National Center for Atmospheric Research (NCAR), which is a major facility sponsored by the NSF under Cooperative Agreement 1852977. E. A. Maroon acknowledges support from the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin–Madison with funding from the Wisconsin Alumni Research Foundation. Keith Lindsay at NCAR aided in ocean model initialization and configuration. Jessica Stevens assisted with data analysis for an early version of this manuscript.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

The Community Earth System Model is freely available from the NCAR but requires registration at www.cesm.ucar.edu/models/cesm1.2. Representative ocean model output needed to reproduce the main figures are publicly available on Open Science Framework using https://doi.org/10.17605/OSF.IO/5ZVWE. Additionally, atmospheric model output for the 150 Tg case (Coupe et al., 2019) is available at https://doi.org/10.6084/m9.figshare.7742735.v1. The full model outputs for all simulations are very large and stored on the PetaLibrary at the University of Colorado, which is not available to the public. However, additional data from these runs can be provided upon request by contacting C. S. Harrison, J. Coupe or C. G. Bardeen.