Nuclear Niño response observed in simulations of nuclear war scenarios

The climate impacts of smoke from fires ignited by nuclear war would include global cooling and crop failure. Facing increased reliance on ocean-based food sources, it is critical to understand the physical and biological state of the post-war oceans. Here we use an Earth system model to simulate six nuclear war scenarios. We show that global cooling can generate a large, sustained response in the equatorial Pacific, resembling an El Niño but persisting for up to seven years. The El Niño following nuclear war, or Nuclear Niño, would be characterized by westerly trade wind anomalies and a shutdown of equatorial Pacific upwelling, caused primarily by cooling of the Maritime Continent and tropical Africa. Reduced incident sunlight and ocean circulation changes would cause a 40% reduction in equatorial Pacific phytoplankton productivity. These results indicate nuclear war could trigger extreme climate change and compromise food security beyond the impacts of crop failure. A nuclear conflict could lead to global cooling followed by a large and sustained El Niño-like response in the tropical Pacific region, according to an evaluation of six scenarios for nuclear conflict with an Earth system model.

T he threat to the climate posed by nuclear war was initially explored during the 1980s nuclear arms race between the United States and the Soviet Union [1][2][3][4][5][6][7] . Contemporary global tensions and existing nuclear arsenals motivate a reexamination of this threat. The detonation of nuclear weapons would cause cities to burn in mass fires, producing soot (black carbon) that would be lofted into the stratosphere, where it would linger for years 7,8 . Like large volcanic eruptions 9 , asteroid impact events 10 , and some proposed geoengineering schemes 11 , the injection of soot aerosols into the stratosphere following nuclear conflict would cause surface cooling and stratospheric heating. Subsequent simulations of nuclear conflict with increasingly complex climate models have confirmed that sudden cooling would follow large injections of smoke into the upper troposphere and lower stratosphere 8,[12][13][14][15][16][17] . Even a smaller regional nuclear war could threaten human society by imposing extreme perturbations to the climate 16 and compromising regional and global food security by harming agriculture 18,19 . Survivors of the initial explosions and firestorms would have to contend with severely diminished growing seasons and agricultural output 8,18 , and would likely turn to the oceans as an additional source of food. Currently, protein from fisheries and marine aquaculture provides 17% of the world's annually consumed animal protein, with far higher proportions locally 20 . Evaluating the ocean response to nuclear war is thus crucial to assessing potential perturbations to overall food security and provides one motivation for examining the tropical Pacific climate response.
The El Niño-Southern Oscillation (ENSO) is the largest naturally occurring perturbation to Pacific Ocean circulation and biogeochemistry, alternating between warm El Niño and cold La Niña events with a period of roughly 3-7 years 21 , with profound impacts on marine productivity and fisheries 22,23 . ENSO is known to be sensitive to both natural and anthropogenic climate influences 24 , providing a pathway to understand the physical and biological impacts of aerosol-driven cooling on the ocean; yet the response of ENSO to aerosol cooling from nuclear war has not been studied beyond a brief mention of El Niño-like tendencies in nuclear war simulations 17 . The climatic perturbation induced by a nuclear war is similar to the effects of a very large tropical volcanic eruption, but it is longer lasting due to the stronger selflofting properties of the soot released into the atmosphere 2,8 . Paleoclimate data [25][26][27] and modeling experiments [28][29][30][31][32][33][34][35][36] have shown that volcanic eruptions increase the probability of an El Niño event in the following years, with a La Niña to follow 27,37 . Likewise, an analysis of geoengineering experiments using stratospheric aerosols to counteract anthropogenic global warming found a slight increase in the probability of an El Niño event in the following year 38 and an El Niño-like shift in the decadal average temperature 39 . The response of the ENSO system to a global perturbation with near-decadal persistence has not been studied; the present nuclear war simulations provide an opportunity to examine this coupled ocean-atmosphere response at its most extreme.
Using a state-of-the-art Earth system model and six different nuclear war scenarios, we demonstrate that global cooling generated by nuclear conflict initiates a large, sustained circulation response in the tropical Pacific. This "Nuclear Niño" exhibits qualitatively similar behavior to El Niño and its response to volcanic aerosol forcing, but is of a much larger magnitude and duration. The strong, almost decadal circulation response that occurs in the largest nuclear war scenario is unprecedented in the modern or paleoclimate ENSO record, highlighting the extreme nature of climate change after a nuclear conflict. While previous work has reported on some aspects of these simulations 8,40,41 , neither the ENSO response nor the associated physical mechanisms have been explored. However, several mechanisms have been put forth to explain the link between volcanic aerosol-driven cooling and El Niño-like conditions. One theory suggests the equatorward shift of the Intertropical Convergence Zone (ITCZ) from the differential cooling of the Northern and Southern Hemisphere 42,43 could weaken the trade winds along the equator and trigger El Niño initiation 32,33,44 . Others have suggested the response could be a consequence of: cooling the Maritime Continent which reduces the zonal Pacific SST gradient 30 , eastwardpropagating Kelvin waves initiated by the cooling of Africa 34 , or the "ocean dynamical thermostat mechanism" where cooling in the eastern Pacific is moderated by changes in the vertical structure of the ocean 27,45 . We conduct sensitivity tests to understand and quantify the contribution of different mechanisms toward the extreme tropical Pacific circulation response in our simulations.
Strong El Niño events are associated with warming sea surface temperatures (SSTs) and reduced upwelling along the equatorial and eastern Pacific, which has historically limited fisheries production in the Peru-Chile and California upwelling systems through nutrient limitation 22,23,46 . The large variability in equatorial Pacific primary productivity is often responsible for driving changes to globally averaged productivity 22 . For example, the 1997-1998 El Niño event caused global primary productivity to drop by 3 Pg of carbon per year, mostly due to reduced productivity in the eastern Pacific, which contributed to the collapse of anchovy fisheries off the coast of Chile 23 . Billions of dollars of economic activity are subject to temperature variability associated with ENSO, including 60% of global commercial tuna catches, which originate in the warm waters of the western and central Pacific and may benefit from warming during El Niño events 23,46 . To help characterize the impact of nuclear war on marine productivity and the associated repercussions for food security we examine changes to net primary productivity (NPP) in the equatorial Pacific. Large reductions in NPP are simulated following nuclear conflict and the mechanisms for this response are disentangled. Despite the strong link between equatorial Pacific NPP and El Niño events, we find that reductions in solar radiation account for much of the biological response.
Results nuclear war scenario simulations. We use the Community Earth System Model (CESM, version 1.3) with the Whole Atmosphere Community Climate Model (WACCM, version 4) as the atmospheric component 47 , hereafter CESM-WACCM4, to simulate the climate response to smoke from nuclear war. The ocean model component used is the Parallel Ocean Program (POP, version 2), which is coupled to the Biogeochemical Elemental Cycling model. This is the same model configuration used by Coupe et al. 8 , which is able to replicate the periodicity of ENSO in terms of historical Niño3.4 region (5°N-5°S, 170°W-120°W) SST spectral behavior, although its amplitude and duration is somewhat too large 47 .
CESM-WACCM4 has also demonstrated an ability to simulate the ENSO response to volcanic eruptions and is a suitable tool for these simulations 48 . Additional model component technical details are provided in the "Methods" section. We focus on a scenario of nuclear war between the United States and Russia as described by Robock et al. 12 and Coupe et al. 8 , simulated as a 150 Tg soot injection in a layer between 300 hPa and 150 hPa of the atmosphere over those two countries over seven days, starting on 15 May. Hereafter, this case is referred to as NW-150Tg. Five additional nuclear war scenarios 40 , with injections of 5-46.8 Tg of soot over India and Pakistan are also used to assess the sensitivity of the response to different radiative forcings. Table 1 shows the number of weapons and average weapon yields used in each scenario.
Nuclear Niño response. Soot aerosols spread globally after they are injected into the upper atmosphere, blocking shortwave radiation and causing global cooling. In NW-150Tg, a nearly 120 W m −2 reduction in global mean monthly downwelling solar radiation is simulated (Fig. 1a), accompanied by a nearly 10°C reduction in global mean monthly surface temperatures over the course of the next two years (Fig. 1b). The spatial patterns of temperature (Fig. 1c) and precipitation (Fig. 1d) anomalies are characterized by extreme cooling and reduced precipitation over most of the globe. In contrast, the mild cooling conditions and increase in precipitation along the equatorial Pacific are indicative of the Nuclear Niño response. The Nuclear Niño response is further illustrated for all scenarios in Fig. 2 and is characterized by a sudden increase in westerly surface wind stress along the western and central equatorial Pacific in the month after the injection (Fig. 2a). This westerly anomaly, favorable for El Niño development and a fundamental feature of the Nuclear Niño response, continues to increase until it peaks during  July-August-September. Westerly wind stress persists for up to seven years in NW-150Tg and the direction of the trade winds is reversed intermittently for three years. In response to anomalously westerly wind forcing, SSTs in the Niño3.4 region begin warming in all simulations despite the cooling effect of the stratospheric aerosol layer ( Supplementary Fig. 1). The warming of the Niño3.4 region occurs after the westerly wind stress anomaly advects warm equatorial waters eastward, deepening the thermocline in the central and eastern equatorial Pacific for several years ( Supplementary Fig. 2). The importance of zonal advection to Nuclear Niño development is apparent from the equatorial Pacific mixed-layer heat budget for the period 6-12 months postinjection ( Supplementary Fig. 3). Within a year, aerosol-driven cooling masks equatorial Pacific warming. Regardless of the actual SSTs, the zonal SST gradient across the Pacific assists in driving the atmospheric response. To filter out large-scale cooling, we subtract the average SST anomalies between 20°S and 20°N from the SST anomalies in the Niño3.4 region (5°N-5°S, 170°W-120°W) to get the relative Niño3.4 SSTs 34 . We find anomalously warm relative SSTs along the equatorial Pacific for seven years in NW-150Tg (Fig. 2b). The relative Niño3.4 SSTs in NW-150Tg, even when normalized to account for higher amplitude El Niño events in our model, exceed the SST response for both the peak of the 1997-1998 and 2015-2016 El Niño events. The strong physical ocean response is exemplified by the change in the zonal sea surface height gradient between the east and west Pacific Ocean, which is typically negative but shifts positive for all of the simulations except the 5 Tg soot case (Fig. 2c). The normal circulation pattern of the equatorial Pacific Ocean is completely disrupted, as upwelling is shut down and in some cases reversed completely (Fig. 2d). In NW-150Tg, the Southern Oscillation Index (SOI, see Supplementary Methods) drops seven standard deviations within 6 months. The rate of the collapse of the Pacific Walker Circulation (inferred by changes in the SOI) in all simulations with more than 16 Tg of soot is unmatched by any event in the historical record, control run or 1200 years of the CESM Large Ensemble pre-industrial run. This indicates the Nuclear Niño response is well outside the bounds of what might be expected due to unforced ENSO variability in the historical record and in the model climatology, even though CESM-WACCM4 simulates stronger ENSO variability than observations. The tendency for El Niño-like conditions is consistent with similar work modeling volcanic eruptions but clearly is more extreme [28][29][30][31][32][33][34][35][36] .
The initial ENSO state can also influence the evolution of ENSO after volcanic eruptions. In the control run from which the initial conditions for all nuclear scenarios were obtained, neutral to La Niña-like conditions were simulated in the first year when the Nuclear Niño is initiated in the nuclear war scenarios. To examine the influence of initial conditions, we used two additional ensembles of the 5 Tg soot case: one starting with neutral conditions and one with a developing El Niño. Regardless of initial conditions, in all three ensemble members an El Niño developed that continued through December-January-February at least. In the case with the developing El Niño, adding 5 Tg of soot enhanced the strength of the response by 25% and the duration by two months compared to the control run. We conclude that the initial ENSO state affects but does not prevent the development of an El Niño in even the smallest war scenario.
An extreme Nuclear Niño response is not simply confined to the largest forcing scenario. The 27.3 Tg soot case with only 18% of the initial soot achieves more than 60% of the wind response of NW-150Tg. The strength and duration of the wind response is correlated with increased aerosol mass, although a saturation point occurs after 27.3 Tg ( Supplementary Fig. 4). During the peak of the event, the 27.3 Tg, 37 Tg, and 46.8 Tg soot cases all have a similar surface wind stress forcing, although the larger cases have a longer duration of El Niño-like conditions. While there is no long-duration El Niño in the 5 Tg case, the ensemble mean zonal wind stress forcing remains anomalously westerly for nearly two years. Reduced equatorial upwelling also contributes to the prolonged circulation response by maintaining a weakened zonal SST gradient in the Pacific for nearly seven years in NW-150Tg, evidenced by the large positive vertical advection term in the ocean mixed-layer heat budget for the Niño3.4 region ( Supplementary Fig. 3). The extended duration El Niño is also related to a lack of deep convection in the western Pacific, allowing anomalous sinking air over the region for many years in the larger forcing scenarios.
Mechanisms contributing to Nuclear Niño initiation. To understand the mechanisms responsible, we focus on the initiation of El Niño conditions in NW-150Tg, which has the largest perturbation and represents what can be described as the potential upper bound of an El Niño event. A series of changes occur in this simulation that are consistent with mechanisms previously invoked to explain El Niño-like conditions following volcanic eruptions, such as differential cooling of the Maritime Continent relative to the central Pacific, strong cooling of the African continent, and an equatorward shift in the latitude of the ITCZ (Supplementary Figs. 5, 6). However, before the westerly trade wind anomaly spreads across the equatorial Pacific, there is an initial La Niña-like response in the eastern equatorial Pacific characterized by strong cooling, an increase in coastal upwelling, and anomalous easterly trade winds ( Supplementary Fig. 5). The faster cooling of South America in contrast to the surrounding ocean is the most likely culprit for this easterly wind anomaly, which is confined to a localized coastal area in the lower troposphere and persists for less than six months. La Niña-like features can precondition the ocean for El Niño development through ENSO recharge dynamics 29,44 , but this short-lived feature cannot account for the multi-year persistence of the Nuclear Niño response. Rather, the dominant driver in NW-150Tg appears to be the suppression of convection over the Maritime Continent region in response to cooling ( Supplementary Fig. 5). This cooling also creates a localized weakening of the zonal temperature gradient in the western Pacific, contributing to the initial relaxation of the trade winds.
To quantify the contribution of different mechanisms toward the Nuclear Niño response, we conduct a series of additional sensitivity experiments using CESM-WACCM4. Following a similar approach as Khodri et al. 34 , cooling is simulated by increasing the surface albedo to reflect the amount of solar radiation equivalent to that intercepted by soot aerosols in NW-150Tg. Three main experiments were performed, where cooling was applied over the Maritime Continent (EXP-MC), tropical Africa (EXP-TA), and the Northern Hemisphere (EXP-NH), with exact borders defined in Table 2. Surface cooling in EXP-MC and EXP-TA is designed to mimic diminished convection over those regions in NW-150Tg which has been suggested to be important for triggering an El Niño response 30,34 . EXP-NH simulates the differential cooling of the Northern and Southern Hemisphere to mimic the equatorward shift of the ITCZ. This calculation is discussed in further technical detail in the "Methods" section. In these experiments, cooling the Maritime Continent (EXP-MC) and tropical Africa (EXP-TA) produced the strongest El Niño-like responses (Fig. 3). In NW-150Tg, the peak in Pacific Ocean zonal surface wind stress anomalies (130°E-120°W, 5°S-5°N) favorable for El Niño development occurs during July-August-September of the first year. Wind stress anomalies for that same period in EXP-MC, EXP-TA, and EXP-NH are 71%, 31%, and 29% of the NW-150Tg forcing, respectively (Fig. 3a). The peak in wind stress forcing in EXP-MC and EXP-TA occurs months later, during November-December-January. During this peak, wind stress in EXP-MC is actually greater than the peak in NW-150Tg. Maritime Continent region cooling is thus the most effective at generating the initial wind response that triggers the peak Nuclear Niño conditions in NW-150Tg. Reduced precipitation in response to cooling over the Maritime Continent produces anomalous sinking air and disrupts the Pacific Walker Circulation (Fig. 4). The amount of cooling over the Maritime Continent in EXP-MC and NW-150Tg are similar by design but there is a greater reduction in precipitation and greater wind stress in EXP-MC that can be explained by a slightly weaker zonal SST gradient due to the lack of aerosol-driven cooling of the central and eastern Pacific Ocean. In terms of Niño3.4 relative SST anomalies, NW-150Tg's peak in October-November-December of the first year is 0.3°C cooler than EXP-MC's February-March-April peak in the following year. Warmer SSTs peaking later in EXP-MC are expected with a later and stronger wind stress forcing. These results show that suppression of convection in this region, which is tightly linked to variations in the Pacific Walker Circulation, is capable of generating wind stress and SST anomalies sufficient to explain the Nuclear Niño in NW-150Tg. In EXP-TA, a strong El Niño-like response is also simulated, supporting the findings of Khodri et al. 34 that the cooling of the African continent contributes to the development of El Niño-like conditions following aerosol-driven cooling. The collapse of the West African monsoon (Fig. 4c) produces a mid-tropospheric heating anomaly, initiating an eastward propagating Kelvin wave along the equator that suppresses convection in its wake. After passing through the Maritime Continent and western Pacific, the waves generate westerly zonal wind stress anomalies. Surface wind stress anomalies in EXP-TA peak at the same time as in EXP-MC, but the peak is only 44% of the magnitude of EXP-MC. Niño3.4 relative SST anomalies are also cooler in EXP-TA  Table 2 for description of the experiments.
compared to EXP-MC, suggesting the Maritime Continent cooling is the more dominant mechanism in this model configuration. However, the spatial extent of eastern Pacific warming in EXP-TA is more impressive than in EXP-MC (Fig. 3e). The cooling of the equatorial Atlantic Ocean in EXP-TA due to advection of cold air from the African continent leads to an Atlantic Niña-like pattern, which has been linked to warming in the equatorial Pacific 34 . Reduced convection over the equatorial Atlantic (Fig. 4c) disrupts the Atlantic-Pacific Walker Circulation, contributing to the additional warming in the eastern Pacific by reducing the trade winds along the coast of western South America more so than in EXP-MC and EXP-NH. While equatorial Atlantic cooling is also present in NW-150Tg, a coastal warming effect is likely to be masked by massive cooling of the South American continent. The smaller magnitude of wind stress and relative SST anomalies in the traditional Niño3.4 region in EXP-TA (Fig. 3) suggests that this mechanism plays a smaller role in initiating the Nuclear Niño than Maritime Continent cooling. Still, the large wind response and impressive warming pattern in EXP-TA is evidence of an important link between tropical African convection and ENSO.
In EXP-NH we increase the albedo of the entire Northern Hemisphere north of 10°N to simulate differential cooling of the hemispheres while avoiding equatorial regions to estimate the contribution of the equatorward shift of the ITCZ that occurs in NW-150Tg ( Supplementary Fig. 6). The areas in EXP-MC and EXP-TA were excluded from albedo modifications in EXP-NH; cooling in these regions is thus due solely to advection from other locations. In EXP-NH, El Niño conditions are simulated, but were the weakest of all the simulations in terms of wind stress. Wind stress forcing in the Pacific in EXP-NH peaks during August-September-October, more in agreement with NW-150Tg, but its peak is only 30% of NW-150Tg. This peak occurs as the ITCZ remains near the equator in defiance of its seasonal progression northward. The proximity of the ITCZ to the equator in the Pacific is enough to generate El Niño-like conditions, but not enough to generate the full Nuclear Niño response.
The three mechanisms investigated here do not constitute a linear decomposition of the NW-150Tg response, as evidenced by the fact that the linear addition of the wind stress forcing in EXP-TA and EXP-MC would be larger than the response in NW-150Tg. This is unsurprising, given the degree of coupling inherent in the tropical Pacific climate system. A saturation effect exists for higher westerly wind stress forcing, meaning that though multiple regions cool in NW-150Tg, there are diminishing returns in the overall wind stress response generated. Additionally, sensitivity tests may not perfectly isolate individual mechanisms because of the potential for exaggerated zonal temperature gradients in EXP-TA and EXP-MC, as well as the advection of cold air into other regions that may influence ENSO. NW-150Tg may also include mechanisms that weaken the El Niño-like response that are not present in the sensitivity experiments, like the cooling of South America. While imperfect, the tests allow a reasonable assessment of the relative contribution of major mechanisms to the initial response. Multiple factors likely work together to initiate the El Niño-like conditions after aerosol-driven cooling, but our results indicate that the cooling of the Maritime Continent generates the strongest initial wind stress response, with a significant contribution from cooling of the African continent and a smaller contribution from the equatorward shift of the ITCZ. After initiation, the multi-year persistence of the Nuclear Niño is driven by sustained trade wind forcing from the aforementioned mechanisms which continue to reinforce the Nuclear Niño even after initiation, in addition to internal Bjerknes feedback processes.
Previous work on volcanic eruptions and ENSO has suggested the ocean dynamical thermostat mechanism, where upwelling in the eastern Pacific moderates aerosol-forced cooling compared to the western Pacific, could contribute to the El Niño-like response following stratospheric aerosol injection 27,28,45 . The dynamical thermostat does not seem to initiate the Nuclear Niño, as the zonal Pacific SST gradient and the westerly zonal wind anomaly initially both increase, in contrast to what would be predicted by preferential cooling of the warm pool in the dynamical thermostat paradigm. In fact, the eastern equatorial Pacific initially cools off more than the rest of the basin ( Supplementary  Fig. 5). Additionally, the vertical thermal structure of the ocean is strongly modified by the massive aerosol-induced cooling (Supplementary Figure 3), preventing the mean vertical temperature gradient from acting as a "thermostat." The physical response of the tropical Pacific to nuclear conflict is clearly extreme, and the long persistence of these circulation changes suggests the potential for sustained and severe impacts on marine productivity and associated food security implications. We next investigate these linkages by diagnosing the post-conflict biogeochemical changes in our simulations.
Ocean biogeochemical responses: NPP reductions. Marine NPP decreases substantially in NW-150Tg (Fig. 5), with potentially severe consequences for regional marine ecosystems. NPP is the net fixation of carbon by photosynthetic phytoplankton, typically expressed as a rate (here, g C m −2 yr −1 ), and drives marine ecosystem and fisheries production [49][50][51] . NPP is equal to the product of phytoplankton biomass and cell division rates, which the latter reflect the combined effects of light, nutrients, and temperature (see "Methods"). Under normal conditions, simulated mean annual NPP in the Eastern Equatorial Pacific biome and the Peru Exclusive Economic Zone (EEZ; Fig. 5b) is high, yielding values of~245 and~280 g C m −2 yr −1 (Fig. 5a), similar to estimates from satellite observations 50 . In the 2-5 years after the 150 Tg soot injection, annual mean NPP declines by 37% to 155 g C m −2 yr −1 and 34% to 186 g C m −2 yr −1 , respectively in the two regions (Fig. 5b). This reduction is even more dramatic on the seasonal scale (Fig. 5d), particularly in the Peru EEZ. Reductions in monthly averaged NPP in the Peru EEZ routinely exceed 50%, and reach as high as 68% post-war, a reduction of 200 g C m −2 yr −1 . The NW-150Tg NPP reductions are many times larger than the natural variability, as are the reductions for smaller nuclear war scenarios (Supplementary Fig. 7). These results indicate that nuclear conflict would stress ecosystems not adapted to events of this magnitude, likely leading to sustained decrease in regional marine productivity, in turn reducing fish catch 23,41,46,49 .
Ocean biogeochemical responses: NPP mechanisms. Reductions in NPP in the equatorial and eastern Pacific are typically observed during El Niño events, but these are much smaller in magnitude and are driven by deepening of the thermocline and reduced upwelling of nutrient-rich water 52,53 , a different mechanism than the dominant one simulated here. In the aftermath of nuclear war, large regional NPP reductions in the Eastern Equatorial Pacific biome and Peru EEZ are primarily driven by a reduction in incident light reaching the ocean's surface, i.e., surface photosynthetically active radiation (PAR). In the Eastern Equatorial Pacific biome ( Supplementary Fig. 8), for example, simulated phytoplankton division rates are co-limited by nutrient, temperature, and light limitation, all of which are further modified by deeper mixed layer depths as a result of surface cooling (Supplementary Fig. 8a) and reduced surface PAR ( Supplementary  Fig. 8b) following the conflict. From a nutrient perspective, equatorial Pacific NPP is only affected by changes to the most limiting nutrient, which is consistently iron throughout all control and nuclear war scenarios. Despite a reduction in equatorial upwelling (Fig. 2d) of nutrient-rich water, iron-dominated nutrient limitation, is actually relieved ( Supplementary Figure 8c) due to the combined effects of deeper mixed layer depths ( Supplementary Fig. 8a), which mix nutrient-rich water up toward the surface, and reduced NPP, which decreases iron uptake. The net result of relieved nutrient stress, however, cannot explain the reduction in NPP. Instead, the simulated reduction in NPP and the associated photosynthetic division rates (Supplementary Fig. 8d) must be driven by changes to temperature and light limitation, which are both exacerbated ( Supplementary  Fig. 8c). The most pronounced effect is on light limitation, which is~50% more limiting in the year following the conflict relative to the control (see Biogeochemical Elemental Cycling model in "Methods" for details on this calculation). Changes to light limitation are dominated by a reduction in PAR across the upper water column, rather than deep mixing causing effective light limitation. Together, the combined effect of all three modified limitation terms is a 25% reduction in population specific division rates averaged over the five years following the conflict (Supplementary Fig. 8d), driving a 15% reduction in biomass (Supplementary Fig. 8e) and 34% reduction in NPP (Fig. 4d) over the same period. When forced with smaller soot injections, the biogeochemical mechanisms leading to equatorial Pacific NPP reductions are qualitatively consistent, but with roughly proportionally reduced magnitude (see Supplementary Fig. 8).

Discussion
The climate response to large-scale cooling suggests a distinct, extreme mode of ENSO which was not previously anticipated. An initial El Niño-like response was expected, as existing research links volcanic eruptions and ENSO, but this type of long duration response has not previously been found in climate model simulations with stratospheric volcanic sulfate aerosols, which have a shorter residence time. The peak radiative perturbation at the surface as a result of the soot aerosols in even the smallest 5 Tg scenario simulation (−10 W m −2 , Fig. 1a) is greater than the 1991 Pinatubo eruption (−6.5 ± 2.7 W m −2 ), but less than the perturbation from the highly impactful 1815 Tambora eruption (−17.2 ± 4.9 W m −2 ) 54 . The long residence time of the soot aerosols, more than three times that of volcanic aerosols across all of the cases, is a key factor in this response, especially for cases with more than 5 Tg of soot.
For years, temperature and precipitation anomaly patterns in the Pacific Ocean are comparable to those observed during strong, basin-wide El Niño events. The atmospheric teleconnections associated with Nuclear Niño most likely differ from a traditional El Niño; however, a "nuclear winter" with plummeting global mean surface temperatures occurs regardless of any influence of the tropical Pacific Ocean. The Nuclear Niño response is terminated following years of aerosol removal, when solar radiation levels rise (Fig. 1a) and the Maritime Continent begins warming. Termination of the response also coincides with the warming of tropical Africa and the return of the ITCZ to its typical oscillatory pattern. Afterward, there is a rebound toward La Niña-like conditions as the western Pacific warms more quickly than the eastern Pacific. This rebound is consistent with modeling work and observations after volcanic eruptions and large El Niño events in general 27,37 , but occurs years later. Ultimately, there is no clear process to mitigate the duration or strength of the Nuclear Niño event until enough aerosols are removed from the stratosphere and large-scale cooling can no longer be maintained.
We have shown that a range of nuclear war scenarios could induce an El Niño-like pattern of unprecedented magnitude, initiated by the cooling of the Maritime Continent, the cooling of tropical Africa, and the equatorward shift of the ITCZ. We provide confirmation that the African land-cooling mechanism brought forward by Khodri et al. 34 likely plays a role in the El Niño-like response to aerosol-driven cooling. However, our results, using a different model and a larger negative radiative perturbation, disagree with other work on the relative importance of African land-cooling compared to other mechanisms. This suggests mechanisms may be dependent on the model and radiative forcing, which merits further investigation with additional climate models.
We have quantified the negative effect of reduced light and surface cooling following an all-out nuclear war on net primary productivity, the base of the marine food web. Smoke-driven reduction in solar radiation drives NPP reduction across the equatorial Pacific following the war. Although these substantial reductions in NPP are not predicated on the El Niño-like circulation response, they are exacerbated by the reduction in upwelling. In many simulations of large volcanic eruptions, an initial short-lived El Niño response also reduced equatorial Pacific upwelling but changes in biological productivity were modeldependent 55 . A robust enhancement of biological productivity was then found in the 3-5 year period after the eruption as La Niña-like conditions favored enhanced upwelling and higher NPP 55 . After nuclear war, the longer residence time and the stronger radiative perturbation of soot aerosols would impose greater stress on marine ecosystems, extending the duration of biologically unproductive conditions in the equatorial Pacific compared to volcanic eruptions. Unlike after volcanic eruptions, our results do not show a prolonged period of enhanced productivity when the El Niño response ends. The stress of lightdriven reductions to NPP will further impact fisheries already responding to changes in SST, oxygen and salinity traditionally associated with El Niño-like climate variability. Ensuing changes in the size, location, and distribution of fish populations could devastate fisheries, particularly in the southeast Pacific Ocean. Beyond the tropical Pacific, a tendency for lower NPP is also simulated, especially in NW-150Tg (see Supplementary Fig. 9).
The implications of significantly reduced ocean productivity, along with previously reported reductions in agriculture (e.g., ref. 19 ), are grave. In the largest India and Pakistan nuclear war scenario used in this study, fatalities are estimated to be near 300 million people 40 , and even if the entire populations of the United States and Russia perished from the direct effects of a nuclear war, more than 7 billion survivors would be left behind across the world. Reduced agricultural output and fish catch following a nuclear war are likely to exacerbate food insecurity globally.
The simulations here add to the body of work showing the potential global consequences of a nuclear war and provide a first look at the mechanisms responsible for potential changes in marine circulation and productivity. Our findings further illustrate the potentially catastrophic global consequences of a nuclear conflict, and highlight the role of coupled climate dynamics in sustaining climate change across multiple years. WACCM4 is a chemistry-climate model with interactive chemistry in the stratosphere, allowing the model to simulate the destruction of ozone that occurs with stratospheric heating. We use the Community Aerosol and Radiation Model for Atmospheres (CARMA), a sectional aerosol model coupled to WACCM4, which treats black carbon aerosols as fractal particles. 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. We examine in total 6 nuclear war scenarios. In the United States and Russia nuclear war scenario, 150 Tg of pure black carbon is injected into the 150 hPa to 300 hPa layer over the United States and Russia as a worst case scenario. The injection begins on May 15th at a rate that decreases linearly over one week. As indicated by Table 1, this injection can occur through a range of conflicts. It is possible with 4000 nuclear weapons with a yield of 100 kilotons or with as few as 3100 weapons with yields up to 500 kt. Other than the 150 Tg US-Russia case, we examine five simulations of India-Pakistan nuclear war scenarios with 5 Tg, 16 Tg, 27.3 Tg, 37 Tg, and 46.8 Tg black carbon injections on May 15th over those two countries. Each simulation is run for 15 years. Three ensembles of the 5 Tg injection case were conducted. Additionally, three control run ensembles were run for 20 years each, which is used as a climatology to compare all variables to typical conditions.

Methods
The India and Pakistan scenarios are based on the work by Toon et al. 40 , which assume different nuclear arsenal sizes for each country within the range of present estimates. Toon et al. 40 assumes both India and Pakistan have 250 nuclear weapons each, based on a projection for 2025, and consulted policy and military experts to determine the most plausible scenarios as possible. It is assumed each weapon hits a unique target in an area important for military or industry. The area of fire ignition for each weapon is based on the area observed in Hiroshima (13 km 2 for a 15 kt nuclear explosion) and is scaled linearly up. The amount of smoke generated as a result is a function of population density, as each person is allocated 11,000 kg of flammable material, consistent with the fuel availability reported in the TTAPS study 2 . The simulations presented here are limited by the nuclear war scenarios and smoke estimates that have been developed to date. The other nuclear states, China, North Korea, the United Kingdom, and France, are excluded in our simulations not because the arsenals possessed by these countries do not pose legitimate threats, but wars involving them would produce amounts of smoke in the range that we are already studying. Future work will include different nuclear war scenarios, but any nuclear war would be catastrophic. We make no claims as to the probability of any nuclear war, and hope that our work keeps it to 0%.
Nuclear Niño mechanism sensitivity experiments. The sensitivity tests to help quantify the contribution of different mechanisms toward the Nuclear Niño response in NW-150Tg were conducted with CESM-WACCM4 by selectively cooling different regions. Following the approach of Khodri et al. 34 , cooling is simulated by increasing the surface albedo to reflect the amount of solar radiation equivalent to that intercepted by soot aerosols in NW-150Tg. To estimate the amount of solar radiation intercepted by soot aerosols in NW-150Tg, the daily shortwave radiation anomaly as a fraction is extracted from NW-150Tg. This is used to calculate the increase in albedo (α) that is necessary to reduce absorbed surface shortwave radiation. We calculate this quantity, dα, for each day and grid point: dα ¼ Àcf *ðSWradiation 150Tg À SWradiation climatologyÞ=SWradiation climatology Swradiation_climatology is the daily climatological value of shortwave radiation using the control run. A correction factor (cf) to eliminate the influence of clouds is calculated by taking the ratio of the monthly net surface solar flux and net clear-sky surface flux. This output was not available at the daily level so it had to be interpolated daily. dα is added to the direct shortwave albedo in CESM within the radiation code. This results in cooling similar to NW-150Tg. Three main experiments were performed, where cooling was applied over the Maritime Continent (EXP-MC), tropical Africa (EXP-TA), and the Northern Hemisphere (EXP-NH), with exact borders defined in Table 2.
Ocean biogeochemistry model. Ocean biogeochemistry in this simulation is treated with the Biogeochemical Elemental Cycling (BEC) model 56 . Global solutions have been widely validated against global data sets and shown to capture basin-scale spatial distributions in primary production, nutrient and chlorophyll concentrations [56][57][58][59][60] in addition to key aspects of oceanic iron 61 and carbon cycling 61,62 . BEC features a single class of zooplankton and three classes phytoplankton: diatoms, small phytoplankton, and diazotrophs. Phytoplankton carbon biomass is resolved independently for each phytoplankton pool and tracked in terms of grid cell concentration. Class specific phytoplankton net population growth is governed by a photosynthetic net primary productivity term and opposed by a loss term. The net primary productivity term is equal to a volumetric specific photosynthetic division rate multiplied by the biomass concentration. This volumetric specific photosynthetic division rate is subject to limitation by temperature, light, and multiple nutrients (N, P, Si, Fe). Each unitless limitation term, computed prognostically for each phytoplankton class, varies from 0-1 and multiplicatively scales the maximum division rate. Lower values are more limiting and translate to reduced growth. Temperature, light, and nutrients are co-limiting. Multi-nutrient limitation is treated following Liebig's Law of the Minimum 63 such that the maximum specific division rate is only scaled by the most limiting (i.e., the smallest) nutrient limitation term.
Surface biomass is not always well correlated with the depth-integrated biomass inventory 64 , particularly when the mixed-layer depth is very deep as seen in this simulation. Accordingly, rather than focusing on only the surface population, we express biogeochemically relevant variables as population averages by computing biomass-weighted, depth averaged value to capture the complete biological signal. These population-averaged variables are robust to variability in community composition, dilution during deep mixing, biomass that may exist below a shallow mixed layer, and the possibility of non-uniform biomass profiles or subsurface maxima. Population averages are averaged over the first 150 m of the water column because storage limitations precluded saving the complete depth profiles for all prognostic tracers. However, population-averaged variables were also computed and analyzed using extrapolated tracer profiles to account for regions with mixed layer depths below 150 m and did not qualitatively change the nature of the results. Finally, to help distinguish between the relative effects of biomass redistribution via deeper mixing and light reduction via reduced incoming surface PAR on light limitation, MLD-only and PAR-only light limitation terms were diagnostically approximated. The MLD-only term was computed by taking biomass profiles from NW-150Tg and putting them in the light limitation field from the Control Run, then computing a biomass weighted, depth averaged light limitation term. The PAR-only term was conversely computed by taking biomass profiles from the control run and putting them in the NW-150Tg light limitation field.

Data availability
Input data and files necessary to reproduce our specific experiments are available at https://doi.org/10.6084/m9.figshare.13469907.v1. Postprocessed model output is available for the 150 Tg simulation at https://doi.org/10.6084/m9.figshare.7742735.v2. Model output for additional simulations can be requested from the corresponding author.
Code availability