Climate change is increasing the risk of a California megaflood.

Despite the recent prevalence of severe drought, California faces a broadly underappreciated risk of severe floods. Here, we investigate the physical characteristics of "plausible worst case scenario" extreme storm sequences capable of giving rise to "megaflood" conditions using a combination of climate model data and high-resolution weather modeling. Using the data from the Community Earth System Model Large Ensemble, we find that climate change has already doubled the likelihood of an event capable of producing catastrophic flooding, but larger future increases are likely due to continued warming. We further find that runoff in the future extreme storm scenario is 200 to 400% greater than historical values in the Sierra Nevada because of increased precipitation rates and decreased snow fraction. These findings have direct implications for flood and emergency management, as well as broader implications for hazard mitigation and climate adaptation activities.

INTRODUCTION

California is a region more accustomed to water scarcity than overabundance in the modern era. Between 2012 and 2021, California experienced two historically severe droughts—at least one of which was likely the most intense in the past millennium (, )—resulting in widespread agricultural, ecological, and wildfire-related impacts (, ) and ongoing drought-focused public policy conversations. Yet, historical and paleoclimate evidence shows that California is also a region subject to episodic pluvials that substantially exceed any in the meteorological instrumental era ()—potentially leading to underestimation of the risks associated with extreme (but infrequent) floods. Observed extreme precipitation and severe subregional flood events during the 20th century—including those in 1969, 1986, and 1997—hint at this latent potential, but despite their substantial societal impacts, none have rivaled (from a geophysical perspective) the benchmark “Great Flood of 1861–1862” (henceforth, GF1862). This event, which was characterized by weeks-long sequences of winter storms, produced widespread catastrophic flooding across virtually all of California’s lowlands—transforming the interior Sacramento and San Joaquin valleys into a temporary but vast inland sea nearly 300 miles in length () and inundating much of the now densely populated coastal plain in present-day Los Angeles and Orange counties (). Recent estimates suggest that floods equal to or greater in magnitude to those in 1862 occur five to seven times per millennium [i.e., a 1.0 to 0.5% annual likelihood or 100- to 200-year recurrence interval (RI)] (, ).

The extraordinary impacts resulting from GF1862 provided motivation for a 2010 California statewide disaster scenario—known as “ARkStorm” (ARkStorm 1.0)—led by the U.S. Geological Survey in conjunction with a large, interdisciplinary team (). The meteorological scenario underpinning the ARkStorm 1.0 exercise involved the synthetic concatenation of two nonconsecutive extreme storm events from the 20th century (). Subsequent analysis suggested that such an event would likely produce widespread, catastrophic flooding and subsequently lead to the displacement of millions of people, the long-term closure of critical transportation corridors (), and ultimately to nearly $1 trillion in overall economic losses (2022 dollars) ().

Meanwhile, a growing body of research suggests that climate change is likely increasing the risk of extreme precipitation events along the Pacific coast of North America (, ), including California (–), and of subsequent severe flood events (, ). The primary physical mechanism responsible for this projected regional intensification of extreme precipitation is an increase in the strength of cool-season atmospheric river (AR) events (–). Previous analyses have suggested that the thermodynamically driven increase in atmospheric water vapor with warming is directly responsible for most of this projected AR intensification [e.g., ()], with the remainder contributed by shifts in regional atmospheric circulation. There is also evidence that increased radiative forcing may result in an eastward shifted expression of atmospheric circulation anomalies associated with both the Madden-Julian Oscillation () and the El Niño–Southern Oscillation (ENSO)–forced component of the Pacific North American pattern ()—both of which would increase the subseasonal variability of cool season precipitation over and near California. Compounding the increase in extreme precipitation associated with AR events are warming temperatures themselves ()—which raise the mean elevation of snow accumulation in mountainous areas (), increase instantaneous runoff rates as rain falls at the expense of snow (), and raise the risk of “rain on snow” events (). Collectively, these previous research findings motivate the question of whether climate change may substantially affect the odds of “low probability but high consequence” flood events.

Here, we describe the overall design and implementation of, as well as results from, “ARkStorm 2.0”—a new severe storm and flood scenario reimagined for the climate change era. Leveraging recent advances in atmospheric modeling by coupling a high-resolution weather model to a climate model large ensemble, we assess the meteorological characteristics of extreme storm sequences (henceforth referred to as “megastorm” events) as well as the subsequent extreme runoff and adverse hydrologic outcomes such meteorological conditions (henceforth, “megaflood” events) would produce under both present-day and warmer future climate regimes. This work builds upon previous research by explicitly considering long-duration (30-day) storm sequences (rather than single-storm events) most relevant to flood hazard management and disaster preparedness, characterizing large-scale ocean and atmosphere conditions associated with such severe storm sequences, and assessing the likelihood of these events over a wide range of potential levels of global warming. We find that climate change has already increased the risk of a GF1862-like megaflood scenario in California, but that future climate warming will likely bring about even sharper risk increases.

RESULTS

Large-scale and regional climate conditions associated with megaflood scenarios

We design two separate megastorm scenarios capable of causing a megaflood in California—one drawn from the recent historical climate (circa 1996–2005; henceforth “ARkHist”) and another from a hypothetical warmer future climate (2071–2080 in the “high warming” RCP8.5 emissions scenario; henceforth “ARkFuture”). Each scenario comprises a multiweek sequence of consecutive severe winter storm events similar to what is reported to have occurred during the peak of the GF1862 event. Specific events are selected by ranking the 30-day cumulative precipitation on a California statewide basis simulated by the 40-member Community Earth System Model Large Ensemble (CESM1-LENS) and subsequently choosing from among the top 3 ranked events in each climate era to dynamically downscale using a high-resolution weather model [the Weather Research and Forecasting (WRF) model v4.3]. Further details can be found in Materials and Methods.

We find that both ARkHist and ARkFuture events occur during simulated warm-phase ENSO (El Niño) years, although the El Niño event that co-occurs with ARkFuture is much stronger [Niño 3.4 sea surface temperature (SST) anomaly = +1.48 K] than that with ARkHist (+0.56 K). Both events have maximum SST anomalies located in the tropical central Pacific (Fig. 1, A and B), which would be consistent with so-called “central Pacific” or “Modoki” El Niño (). Warm (positive) SST anomalies are also present in the western Bering Sea and Sea of Okhotsk, as well as along the immediate California coast, in both cases. In addition, a broad region of negative sea level pressure (SLP) anomalies is centered over the Gulf of Alaska and adjacent portions of western North America—consistent with traditional El Niño teleconnections—although the zone of negative SLP anomalies extends farther westward across the North Pacific in ARkHist.

Fig. 1.

Large-scale conditions during California megastorm scenarios.

(A and B) Mean SST anomalies (color contours, K) and mean SLP (hPa) anomalies (dashed/solid contours) during ARkHist (A) and ARkFuture (B). SST and SLP are detrended before anomaly calculation using monthly data from each corresponding CESM1-LENS member (baseline period 1980 to 2005 for ARkHist; 2060–2090 for ARkFuture); solid (dashed) SLP contours denote positive (negative) anomalies in increments of 2 hPa. (C and D) Composite instantaneous vertically IVT (kg m−1 s−1) for all hours in which California mean precipitation exceeds 1.5-mm ARkHist (C) and ARkFuture (D) using WRF 81-km simulations. Mean 30-day 500-hPa geopotential height (GPH, detrended) anomalies (color contours, m) and mean absolute 850-hPa wind vectors (m/s) (black arrows) during ARkHist (E) and ARkFuture (F).

We acknowledge, however, that these large-scale patterns and associations with ENSO are drawn from only two individual scenario instances, and we cannot determine from this analysis alone whether these relationships are robust across a wider range of potential megastorm events. To offer a more systematic assessment, we consider the top 4 ranked 30-day California precipitation events in the CESM1-LENS historical and warmer future snapshot periods (fig. S1). We find that all eight such events are associated with anomalously warm conditions in the tropical Pacific Ocean, and Niño 3.4 SST anomalies are uniformly positive (+0.33, +0.56, +2.28, and + 1.56 K for the top 4 historical events and + 1.17, +1.95, +1.48, and + 1.39 K for future events, respectively, using detrended SST). However, it has recently been demonstrated that dynamic ENSO indices can better capture the spatial diversity of ENSO events and their subsequent western U.S. hydroclimate teleconnections (). We thus calculate the ENSO Longitude Index (ELI)—an ENSO metric that tracks the average longitudinal position of ENSO-associated deep convection and accounts for the nonlinear response of convective activity to SST (). As with Nino 3.4 SST anomalies, all eight such events are again associated with anomalously warm conditions in the tropical Pacific Ocean, but ELI values more clearly illustrate a wider range of ENSO spatial variability and dynamical intensity (ELI = 169.9°E, 171.6°E, 185.1°E, and 181.5°E for the top 4 historical events and ELI = 174.2°E, 181.0°E, 176.8°E, and 179.1°E for future events, respectively, using detrended SST).

Using the ELI categorizations defined in (), this suggests that two of four events each in the historical and future simulations occur under “strong El Niño” conditions (ELI ≥ 179°E), and one of four historical and two of four future events occur under “moderate El Niño” conditions (170°E ≤ ELI < 179°E), with the final historical event falling nominally under the “moderate” threshold. Collectively, seven of eight historical and future potential California megastorm events occur under moderate or strong El Niño conditions as defined by the ELI (eight of eight, if rounding to the nearest degree of longitude). These findings strongly suggest that there is a substantially elevated likelihood of month-long storm sequences capable of producing very large precipitation accumulations during moderate to strong El Niño conditions and that the conspicuous anomalous deepening of the Gulf of Alaska low present in most of these eight events (fig. S3) is plausibly linked to El Niño teleconnections [which would be consistent with ()].

Much previous work has focused on the critical role AR storms (“ARs”) play in California hydroclimate—both as beneficial bolsterers of water supply and as the cause of hazardous floods (–). Composite analysis of 30-day averaged vertically integrated water vapor transport (IVT) and animations of IVT over the 30-day scenarios (movies S1 and S2) confirm that ARs are the primary storm mode during both ARkHist and ARkFuture (Fig. 1, C and D) scenarios, with a well-defined moisture transport axis extending northeastward from just north of the Hawaiian Islands to central California. This alignment is suggestive of 30-day mean storm trajectories capable of entraining large quantities of subtropical moisture (i.e., a “Pineapple Express”–type pattern), although with considerable upstream longitudinal extension of the IVT corridor westward of Hawaii (particularly in the future scenario; Fig. 1, C and D). This overall zonal pattern (but with localized meridional flow near California) is consistent with that recently associated with “AR families” occurring during El Niño conditions (), which tend to be characterized by a strengthened subtropical Pacific jet stream and a persistently anomalous Gulf of Alaska cyclone that together favor long-duration periods of successive AR activity across California. While the general spatial structure of IVT is similar for both scenarios, ARkFuture exhibits mean 30-day composite IVT values that are ~25% higher than ARkHist.

Both severe storm sequences are associated with strong westerly (zonal) winds throughout nearly the entire atmospheric column (fig. S4), with a pronounced vertical maximum of ~60 m/s located around jet stream level (200 to 250 hPa) between 30°N and 35°N. Zonal winds are stronger in ARkFuture, especially in the upper troposphere (by >10 m/s above ~400 hPa). Analysis of 500-hPa geopotential height fields (Fig. 1, E and F) indicates that both events are associated with a broad region of negative mid-tropospheric height anomalies over the North Pacific to the west of California, although the negative height anomaly is more localized to the northeastern Pacific in ARkFuture. This suggests that both ARkStorm scenarios are associated with a robust Pacific jet, which is dynamically consistent with the eastward extension of the wintertime Pacific jet associated with both El Niño (Fig. 1, A and B) [e.g., ()] and climate change [e.g., ()], although the 30-day mean low-level (850-hPa) flow pattern exhibits a slightly more zonal pattern (with less of a meridional component over the northeastern Pacific) in ARkFuture relative to ARkHist. Visual inspection of movies S1 and S2 further confirm that both 30-day scenario storm sequences are characterized by the occurrence of multiple deep extratropical cyclones just west of or over California, which is consistent with recent results in (), which found that AR-associated precipitation in the San Francisco Bay Area increased more for ARs directly associated with extratropical cyclones than those without.

We also find that composite atmospheric instability is relatively high during both ARkStorm scenarios. A 30-day composite convective available potential energy (CAPE) exhibits a broad region of >300 J/kg west of the northern California coast during ARkHist, with an even wider region of CAPE (>300 J/kg) (and locally >400 J/kg) in ARkFuture (fig. S5). The values might be unremarkable in a different geographic context, but in coastal California, ARs are typically associated with primarily stratiform or dynamically forced precipitation, and California ARs tend to be characterized by moist-neutral (versus conditional unstable) vertical profiles (). Modest increases in atmospheric instability have been associated with outsized impacts during certain historical California storm events, increasing the risk of flash flooding/debris flows () and severe wind gusts () (fig. S6).

Cumulative and extreme precipitation

In both ARkHist and ARkFuture, 30-day cumulative precipitation is extremely high. In ARkHist, we find broad regions exceeding 500 mm of cumulative precipitation, with widespread areas exceeding 1000 mm in the Sierra Nevada (SN) and more isolated pockets exceeding 1000 mm in the Coast Ranges, Transverse Ranges, and far southern end of the Cascade Range (domain maximum of ~2150 mm; Fig. 2A). In ARkFuture, spatial patterns of event total precipitation are similar but are uniformly characterized by heavier accumulations, with broad areas in both northern and southern California exceeding 700 mm and widespread areas in the abovementioned mountain areas above 1400 mm (domain maximum of ~3200 mm; Fig. 2B). We note that these values are comparable to maximum precipitation informally reported during the GF1862, which exceeded 2500 mm in at least two locations on the SN western slope over a slightly longer (~40-day) period (). In general, cumulative precipitation in ARkFuture is between 35 and 60% higher than in ARkHist for northern and central California (although locally >80% higher), with lesser increases in far southern California (fig. S7, A and B). On a statewide average basis, 30-day precipitation is ~45% higher in ARkFuture.

Fig. 2.

Precipitation associated with California megastorm scenarios.

(A and B) Cumulative 30-day precipitation (mm) during ARkHist (A) and ARkFuture (B). (C and D) Cumulative number of heavy precipitation days (days with precipitation > 20 mm/day) during ARkHist (C) and ARkFuture (D). (E and F) Cumulative number of heavy precipitation hours (hours with precipitation >10 mm/hour) during ARkHist (E) and ARkFuture (F). (G and H) Time series depicting hourly precipitation (mm/hour) on a cumulative California statewide basis during ARkHist (G) and ARkFuture (H). Data depicted in all panels are from the innermost 3-km WRF domain.

Although absolute increases in cumulative precipitation are highest in mountainous areas (fig. S7A), relative increases in event total precipitation are greatest in areas that are not prone to orographic enhancement of precipitation during prevailing southwesterly winds (fig. S7B). Thus, some of the largest relative increases in precipitation (locally >80%) instead occur in regions that are less historically accustomed to receiving extreme precipitation during these events, such as inland valleys and otherwise wind-shadowed areas, which is consistent with earlier work ().

Both ARkStorm scenarios are also notable for their very high precipitation intensities. We quantify this on several time scales, focusing on the frequency (over the 30-day scenario periods) with which precipitation intensity exceeds fixed daily and hourly thresholds [the number of days with precipitation > 20 mm/day and the number of hours with precipitation > 10 mm/hour, henceforth “heavy precipitation days” (HPDs) and “heavy precipitation hours” (HPHs)]. In ARkHist, we find that nearly all coastal areas experience at least 8 (of 30) days with precipitation exceeding 20 mm, and most mountain areas exceed 14 such days (except the Transverse Ranges in southern California, Fig. 2C). In ARkFuture, we find a sharp increase in the number of HPDs, especially in northern and central California, where most coastal areas exceed 16 (of 30) HPDs and most mountain areas exceed 20 such days (Fig. 2D and fig. S7, C and D). In some small pockets in the northern SN and far southern Cascades, all 30 days of the ARkFuture scenario are HPDs. HPD increases are substantially smaller in magnitude across southern California (mostly on the order of one to five additional days) but still nearly ubiquitous (fig. S7C).

Because of their particular relevance in the context of flash flood and debris flow risk (), we specifically consider the occurrence of short-duration precipitation extremes in both ARkStorm scenarios. We find that the highest number of such hours occur in orographically favored areas, with the highest frequency of occurrence in the southern California Transverse Ranges and the Feather River watershed in the northern SN during ARkHist (Fig. 2, E and F). In ARkFuture, we report large and widespread increase in the occurrence of HPHs across essentially the entire domain. The largest increases [+25 to 40 cumulative hours (fig. S7, E and F)] occur broadly across the SN and (locally) in Santa Lucia Mountains—shifting the domain-wide maximum in HPH from southern to northern California. We find large relative increases (~200 to 300%) in the frequency of HPH and a large increase in the spatial extent of affected regions in ARkFuture. On a statewide average basis, we find that the frequency of HPH is ~220% higher in ARkFuture versus ARkHist (Fig. 2, G and H). Oakley et al. () conducted a literature review on published hourly rainfall rates in California and/or similar Mediterranean climate regions thought to be sufficient to trigger shallow landslides and debris flows in susceptible terrain, noting a range (5 to 20 mm/hour) that encompasses our HPH threshold (10 mm/hour) in the present study. These findings, therefore, likely have large implications from a flash flood and debris flow risk perspective.

California-wide average cumulative precipitation during the 30-day periods encompassing both extreme storm sequence scenarios represents a considerable fraction of the total annual [October-September water year (WY)] precipitation occurring during both ARkHist (~447 mm or 46% of the WY total) and ARkFuture (~586 mm, of 40% of the WY total). Compared to the climatological mean WY precipitation across all 40 ensemble members during the baseline periods (1996–2005 and 2071–2080, respectively); however, these events represent an even larger fraction of average annual precipitation—60% of WY precipitation in ARkHist and 71% of WY precipitation in ARkFuture. This also means that both the ARkHist and ARkFuture occur during anomalously wet WYs overall (32 and 77% wetter than the contemporaneous averages in ARkHist and ARkFuture, respectively). This would be dynamically consistent, from an ENSO teleconnection perspective, with the strong relationship between moderate to strong El Niño events (as characterized by the ELI) and anomalously wet cool-season conditions in California (). It also has significant implications from a potential flood hazard perspective, as soil conditions are likely to be more saturated than average during anomalously wet WYs, likely amplifying runoff and further elevating the risk of flooding.

To systematically contextualize the precipitation-related results arising from these two specific downscaled extreme storm scenarios drawn from CESM1-LENS relative to all top-ranked 30-day precipitation events in multiple large ensembles—including the CanESM2, GFDL-CM3, and CSIRO-Mk3.6 ensembles [as described in ()]. We conducted an intercomparison of these events during the historical and future study periods. We found that of the top 4 ranked megastorm events (as quantified by California-wide cumulative 30-day precipitation), all 16 events across the four single-model large ensembles have larger cumulative precipitation in the warmer future scenario versus their counterparts drawn from cooler historical climate snapshot period (fig. S1). We further show that hourly precipitation maxima are also higher in future versus historical megastorm events in all four large ensembles (fig. S2).

We also note that there are substantial differences across the large ensembles regarding the absolute magnitude of the 30-day precipitation associated with the top four ranked storm sequences, with CESM1-LENS exhibiting the largest precipitation accumulations (fig. S1). However, a direct comparison between these absolute precipitation values is not possible in this context because of the widely differing number of ensemble members and potential and biases in the representation of extreme precipitation in specific models. Nevertheless, we emphasize that the overall consistency of the response of both 30-day cumulative and hourly precipitation in the warmer future versus cooler historical megastorms, in relative terms within each respective large ensemble, suggests that many of the key conclusions drawn from the two synthetic case studies drawn from CESM1-LENS and emphasized in this analysis are likely to be generalizable.

Precipitation phase, freezing level height, and snow water equivalent

The heaviest precipitation during both ARkStorm scenarios occurs over mountainous terrain—particularly in the SN—and a substantial fraction of that high elevation accumulation falls in the form of snow. In ARkHist, a substantial fraction of the higher elevation portions of the SN receives more than 1000 mm (Fig. 3A) of snow water equivalent (SWE) over the 30-day event (yielding a domain maximum of 7.7 m of accumulated snowfall). Estimates of peak on-the-ground SWE range from around ~300 mm in the southern Sierra to 470 mm in the central Sierra (fig. S8), with even higher maxima over localized mountain peaks (Fig. 3). This extremely heavy snowfall would likely be highly disruptive to infrastructure and emergency response activities.

Fig. 3.

Snowfall and surface runoff associated with California megastorm scenarios.

(A and B) Cumulative 30-day gross SWE (mm) during ARkHist (A) and ARkFuture (B). (C) Difference in cumulative SWE (mm) between ARkFuture and ARkHist. (D and E) Mean snow fraction (snow-to-rain ratio, in percent) during ARkHist (D) and ARkFuture (E). (F) Difference (%) in mean snow fraction between ARkFuture and ARkHist. (G and H) Mean freezing level (m) during ARkHist (G) and ARkFuture (H). (I and J) Time series depicting hourly surface runoff (mm/hour) on a cumulative California statewide basis during ARkHist (I) and ARkFuture (J). Data depicted in all panels are from the innermost 3-km WRF domain.

In ARkFuture, we find that the event-averaged precipitation phase changes from primarily snow to primarily rain at low to mid-elevations (~1200 to 2000 m) but remains primarily snow at very high elevations (≥2500 m) in the SN (Fig. 3, D and E). This results in a spatial dipole pattern of SWE changes, with large (>50%) SWE decreases at lower elevations but large SWE increases at the highest elevations (≥3000 m) of the SN and southern Cascades (locally >50%, yielding cumulative total SWE as high as 1800 mm and a domain maximum of 10.4 m of accumulated snowfall) (Fig. 3, B and C). Further, there is a stark contrast between the large SWE and snow-to-rain ratio decrease in the northern SN versus a substantial SWE increase and lesser snow-to-rain ratio decrease in the southern SN (Fig. 3F) (likely because of lower elevations in the northern Sierra). We report widespread increases in the mean atmospheric freezing level height during ARkFuture (statewide freezing level of ~2230 m for the 30-day window) versus ARkHist (freezing level of ~1940 m; Fig. 3, G and H)—supporting prior studies finding that warmer temperatures during future extreme storm events will fundamentally alter mountain hydrology and subsequent watershed response [e.g., () and ()].

Very large increase in cumulative and peak runoff during ARkFuture

We find that both ARkStorm scenarios are likely to generate very high runoff across a wide range of watersheds and topographies. Projected increases in ARkFuture runoff, however, are widespread and extremely high in magnitude. On a statewide basis, peak runoff during ARkFuture is more than double that during ARkHist (Fig. 3, I and J). In certain key watersheds, however, the relative differences are even larger: In all three SN subregions, the peak runoff is 200 to 400% higher in ARkFuture (fig. S9). A ~100% increase in peak runoff is also observed in the South Coast and Cascade subregions, with a 60% increase along the North Coast.

Event total cumulative runoff increases are similarly large, with increases of 100% or more across most of the SN western slope, the southern Cascades, the Santa Lucias, and also in several major urban areas with a high impervious surface fraction (including the Los Angeles, Sacramento, and San Jose metropolitan areas; fig. S9B). Even greater fractional increases are found for extreme runoff periods (defined as hours with surface runoff of >10 mm/hour; fig. S9, C and D), which increase from being almost negligible in ARkHist (generally three or fewer total hours, except in the Los Angeles Basin) to being widespread across nearly all of California’s major urban areas and mountain ranges (with many locations experiencing >10 such extreme runoff hours). In addition, we find that runoff efficiency during ARkFuture relative to ARkHist (measured as the ratio of total 30-day runoff to 30-day precipitation) increases by ~50% (from ~0.19 to ~0.29)—suggesting that a considerably higher fraction of precipitation is likely to immediately contribute to potential food risk in the warmer future scenario.

Given the geographic concentration of numerous critical pieces of water and flood management infrastructure on the western slopes of the SN Mountains and in California’s Central Valley, we conduct additional analysis focused on the Sacramento and San Joaquin River watersheds that encompass these regions [as defined by their respective U.S. Geological Survey (USGS) Hydrologic Unit Code (HUC) footprints; fig. S10]. We find large and ubiquitous increases in the upper tail of the empirical distribution of both precipitation and surface runoff at both hourly and 24-hour temporal aggregations in ARkFuture relative to ARkHist, although the relative increases are larger for the San Joaquin basin than the Sacramento Basin (Fig. 4). Here, again, we find that the relative increases in the uppermost tail of the surface runoff distributions are much larger than that of the precipitation distributions. At the 24-hour aggregation level, the upper tail of the surface runoff distributions is largely nonoverlapping in both basins (Fig. 4, G and H)—with virtually no overlap at all in the San Joaquin basin during ARkFuture relative to ARkHist. This points to the potential for historically unprecedented surface runoff regimes during future extreme storms in a strong warming scenario—especially in the watersheds draining the western slopes of the central and southern SN, with major implications for operation of critical water infrastructure in these regions.

Fig. 4.

Upper tail of precipitation and surface runoff distribution for Sacramento and San Joaquin River watersheds.

Empirical histograms depicting the cumulative number of hours (over the 30-day scenarios) at or above specific precipitation [purple bars (A to D)] and surface runoff [blue bars (E to H)] thresholds (in units of mm/hour) at two levels of temporal aggregation (1 hour and 24 hours) for two key California watersheds as outlined by HUC Subregion 1802 (the Sacramento River watershed) and HUC subregion 1804 (the San Joaquin River watershed). Data are drawn from the WRF 3-km domain for ARkHist (darker bars) and ARkFuture (lighter bars) and are calculated as average values for each entire watershed. Values less than 1 mm/hour for precipitation and 0.5 mm/hour for surface runoff are excluded from this upper tail analysis.

We attribute these notably high increases in runoff, which greatly exceed fractional increases in precipitation, to the nonlinear hydrologic effects of increasing both total precipitation (via increased AR intensity) and decreasing the snow-to-rain fraction (due to AR warming and the solid-to-liquid phase change of precipitation). This so-called “double whammy effect,” whereby both the volume of precipitation falling on watersheds and the fraction of that precipitation that immediately becomes runoff at higher elevations increases substantially, can be responsible for unexpectedly large increases in runoff volume (). We also suggest that there is arguably a “triple whammy” effect at play in the case of ARkFuture: In addition to the previous two factors, there is evidence for multiple intense “rain on snow” events () in both scenarios (Fig. 3, G and H) that correspond temporally with event-maximum runoff peaks (Fig. 3, I and J). However, we acknowledge that antecedent hydrologic conditions—particularly soil moisture and the extent/moisture content of snowpack leading up to the event—could potentially have large influences on simulated runoff and ultimately on potential flood risks. In this analysis, we only consider the specific antecedent conditions that were actually present in the respective large ensemble members leading up to the simulated events. Although a comprehensive assessment of the various antecedent hydrological contributors to surface runoff is beyond the scope of the present manuscript, more systematic assessments will be conducted in later stages of the ARkStorm 2.0 project.

Megaflood risk increases robustly as function of climate warming

We assess the cumulative and annual likelihood of a 30-day megastorm sequence capable of causing a California megaflood and find that both increase strongly as a function of climate warming. On a high warming emissions trajectory (RCP8.5), we find that the cumulative likelihood of an ARkHist level event begins to accelerate after the year ~2020 period, with corresponding accelerations becoming apparent earlier (~2000) for lesser (50-year RI) and later (~2030) for higher magnitude (200-year RI) events (Fig. 5A).

Fig. 5.

Climate change and California megastorm risk.

(A) Cumulative occurrence of extreme 30-day precipitation accumulations on a California statewide basis as simulated by the CESM1-LENS ensemble. The three blue-green curves denote cumulative occurrence of events equal or greater in magnitude to the ARkHist scenario, as well as for events with approximate RIs of 50 and 200 years. Data are drawn from the historical CESM1-LENS simulations for 1920–2005 and from the RCP8.5 scenario for 2006–2100. (B) Annual likelihood of extreme 30-day cumulative precipitation events as a function of projected global mean surface temperature (GMST; K) anomaly across the 40-member ensemble. Blue-green curves correspond to definitions in (A). GMST anomaly is defined relative to a baseline calculated from the CESM1-LENS preindustrial control run, and both annual likelihood and GMST are smoothed on a 30-year running mean basis.

To accommodate the various Earth system and sociopolitical uncertainties that complicate future predictions of possible greenhouse gas emission trajectories and to facilitate direct comparison with various proposed targets linked to specific planetary warming levels, we conduct further analysis to estimate changes in megastorm risk as a function of the warming itself. We find that the annual likelihood of an ARkHist level event increases rapidly for each 1°C of global warming [by ~0.012/year per degree C from a baseline of ~0.01/year]; Fig. 5B) and that this approximately linear relationship (P < 0.001) appears to hold even at very high levels of warming (~+4°C). We find that climate change to date (as of 2022) has already increased the annual likelihood of an ARkHist event by ~105% relative to 1920 in the CESM1-LENS ensemble and of an even higher magnitude (200-year RI) event by ~234%. This finding is consistent with prior work reporting progressively larger increases in projected extreme precipitation events for increasing event magnitudes [e.g., ()]. We further find that by ~2060, on a high emissions trajectory, the annual likelihood of an ARkHist level event increases by ~374% and by ~683% for a formerly 200-year RI event. These statistics represent notably large increases in risk of California megastorm events due to climate change, as they transform an event that previously would have occurred once every two centuries into one that may occur approximately three times per century.

DISCUSSION

Our analyses suggest that the fundamental characteristics of the plausible worst-case California megafloods of the future will be familiar: Similar to their contemporary and historical counterparts, they will be characterized by a week-long sequences of recurrent, strong to extreme ARs during the cool season and coinciding with a persistently strong Pacific jet stream. Yet, we also find evidence of some critical differences: Future extreme storm sequences will bring more intense moisture transport and more overall precipitation, along with higher freezing levels and decreased snow-to-rain ratios that together yield runoff that is much higher than that during historical events. In addition, we find even larger increases in hourly rainfall rates during individual storm events, which have high potential to increase the severity of geophysical hazards such as flash flooding and debris flows. This is especially true in the vicinity of large or high-intensity wildfire burn areas, which are themselves increasing due to climate change () and yielding large increases in associated compound hazards ().

An extensive body of existing research has linked climate change to increasingly extreme precipitation events [e.g., (–)], even in locations where changes in mean precipitation are nonrobust (, ). There is further evidence that climate warming increases the intensity of ARs in many regions (), including California (, ). The strongest ARs are expected to strengthen considerably at the expense of the weakest—shifting the balance from “primarily beneficial” AR events to “primarily hazardous” ones ()—an intensification brought about primarily via the direct thermodynamic effect of warming ().

Our analysis goes beyond these prior works to demonstrate that climate change is robustly increasing both the frequency and magnitude of extremely severe storm sequences capable of causing megaflood events in California. Our analysis suggests that the present-day (circa 2022) likelihood of historically rare to unprecedented 30-day precipitation accumulations has already increased substantially and that even modest additional increments of global warming will bring about even larger increases in likelihood. Critically, this finding means that existing international emissions policies, which are estimated to yield cumulative warming of well over 2°C (), will entail large further increases in the likelihood of a California megaflood event. We further find that all of the most intense 30-day megastorm events in the CESM1-LENS ensemble occur during moderate to strong ENSO warm phase (El Niño) conditions—both in the historical and warmer future scenarios—suggesting that these events may potentially exhibit some degree of predictability at seasonal scale. For these reasons, we emphasize that recognizing and mitigating the societal risks associated with this subtly but substantially escalating natural hazard is a critically important consideration from a climate adaptation perspective.

Recent evidence suggests that increases in western United States flood risk caused by anthropogenic warming may have been counteracted in recent decades by natural variability, but that further warming and shifts in natural variability will eventually “unmask” this accumulated increase in regional flood risk (). Additional work suggests that the response of flood risk to climate change is likely to exhibit threshold behavior, at least in certain climatological and hydrological regimes (), with a precipitation extremeness threshold dictating whether flood risk decreases (for smaller events, due to the antecedent soil aridification effect of warming temperatures) or increases (for the largest events, due to the overwhelming effect of large increases in precipitation intensity). Both of these considerations are especially germane to California—a region where most contemporary public policy and climate adaptation efforts emphasize drought and wildfire risk due to lack of recent experience with widespread severe floods. Collectively, the findings from previous work and this study illustrate the growing urgency of planning for and mitigating the hazards from potentially catastrophic floods in California in a warming climate.

The extreme storm scenario development and subsequent analyses described here represent the first phase of the broader ARkStorm 2.0 exercise, which is eventually expected to encompass a full suite of follow-on hydrologic and inundation modeling, hazard assessments, and tabletop disaster response exercises. We plan to work with local, regional, and federal stakeholders to integrate quantification of physical hazards resulting from an “ARkStorm”-level event in California within disaster resilience and climate adaptation frameworks. Our initial atmospheric modeling results presented here demonstrate that extremely severe winter storm sequences once thought to be exceptionally rare events are likely to become much more common under essentially all plausible future climate trajectories—suggesting that 20th century hazard mapping, emergency response plans, and even physical infrastructure design standards may already be out of date in a warmer 21st century climate. Still, region-wide and high-resolution runoff inundation modeling capable of accounting for the effects of various active and passive flood management infrastructure will be required to fully quantify the extent of flood-related hazards and associated societal impacts resulting from these two ARkStorm 2.0 scenarios, and these simulations are actively being planned for the project’s future phases.

Yet, potential solutions to increasing flood risk do exist. Examples of climate-aware strategies that have the potential to mitigate harm during a 21st century California megaflood include floodplain restoration and levee setbacks, which would lessen flood risk in urban areas while offering environmental cobenefits (); forecast-informed reservoir operations, which would afford reservoir operators greater flexibility in the face of uncertainty (); and revised emergency evacuation and contingency plans that accommodate the possibility of inundation and transportation disruption extending far beyond that which has occurred in the past century. Some of these interventions—such as flood-managed aquifer recharge—even have the potential to reduce flood damages while simultaneously improving resilience to future regional droughts (). Ultimately, our hope is that the analysis described here can serve as a geographically portable framework for scenario-based emergency response and regional adaptation endeavors in the climate change era, both within and beyond California.

MATERIALS AND METHODS

Overall ARkStorm 2.0 scenario design

ARkStorm 2.0 is a wide-reaching extreme storm and flood scenario for California that seeks to build upon previous disaster contingency and emergency response planning efforts. This endeavor is intended to build upon previous efforts in the original ARkStorm exercise (ARkStorm 1.0), which was completed in 2010 () and involved a broad consortium of local, state, and federal agencies. It was found that the hypothetical storm scenario used in ARkStorm 1.0 would have produced widespread, deep inundation of a large fraction of the Sacramento and San Joaquin valley floors, as well as widespread, life-threatening flooding in other highly populated parts of California. Total economic losses (the sum of direct damages and indirect losses due to business and economic disruption) were projected to exceed $750 billion [2010 dollars ()]. This would be equivalent to approximately $1 trillion in 2022 dollars, making it the most expensive geophysical disaster in global history to date. Partly for this reason, this hypothetical event was informally dubbed California’s “other Big One:” Such a flood event in modern California would likely exceed the damages from a large magnitude earthquake by a considerable margin.

In ARkStorm 1.0, the scenario design involved the artificial concatenation of two of the most intense individual storm sequences in the observed 20th century climate [from January 1969 to February 1986; ()], with additional manual adjustments to the persistence of individual ARs to amplify cumulative precipitation totals. Historical atmospheric reanalysis data were used to obtain boundary conditions for simulating these concatenated events using the Weather Research and Forecasting Model (v3.0.1) at spatial resolution ranging from 2 to 6 km across California. Precipitation and other variables from this single simulation were then used to estimate flood and other related impacts.

In ARkStorm 2.0, we update and upgrade the methods used in ARkStorm 1.0 in several fundamental ways. First, we use a hypothetical extreme event selection method that is both systematic and internally consistent from an atmospheric dynamical perspective: Rather than artificially concatenating multiple historical events, we leverage the large sample size afforded by large ensemble climate model simulations to draw upon a much wider range of physically plausible event sequences that are available by considering the roughly century-long observational record alone (and we make no manual adjustments to storm sequencing). Second, we use a newer and more sophisticated weather model (WRF V4.3) with generally higher spatial resolution (3 km across all of California and adjacent regions). Last and most critically, we design and implement two separate scenarios—ARkHist and ARkFuture—with the combined aim of comparing a “lesser” present era severe storm sequence to a much more intense but physically plausible future sequence amplified by climate change. The overall approach of embedding a high-resolution weather model within existing climate model large ensemble simulations is similar to that described in () and has the dual advantage of not only expanding the statistical sample size of physically plausible but observationally rare or unprecedented precipitation events (in CESM1-LENS) but also attaining the high degree of physical realism afforded by simulating extreme ARs in a high-resolution setting ().

Selection of specific extreme storm sequences

Both ARkHist and ARkFuture are intended to capture multiweek sequences of discrete severe storm events that produce extremely high cumulative precipitation over a 30-day period. The use of a 30-day accumulation period is motivated by the desire to conduct a realistic emergency management contingency exercise as part of ARkStorm 2.0 and the prior knowledge that multiple successive storm events often challenge infrastructure and response systems to a greater degree than shorter-duration events. We first calculate the cumulative 30-day precipitation for the state of California from all 40 ensemble members from the CESM1-LENS () from two decade long “snapshot” intervals during which high-frequency (6 hourly) data are available for dynamical downscaling: 1996–2005 (using the historical scenario, which aims to replicate real-world aerosol and greenhouse gas climate forcings) and 2071–2080 (using the RCP8.5 scenario, which assumes continued rapid growth of greenhouse gas emissions over the 21st century).

Among the available global climate model large ensemble datasets, CESM1-LENS stands out with its comprehensive suite of three-dimensional, high-frequency (6 hourly) atmospheric variables, which provide the forcing conditions required for dynamical downscaling simulations. We note that, while it might otherwise be desirable to sample from a wider time period than the two specific decades included in these snapshots, these are the only two such intervals for which a comprehensive suite of three-dimensional, high-temporal frequency (6 hourly) atmospheric conditions were retained in the original CESM1-LENS experiment, and so, it is not possible to conduct high-resolution WRF simulations during other intervals because of the unavailability of needed initial and boundary conditions. However, as the snapshot periods include data from 40 independent ensemble members initialized decades before the assessment period—each with their own sequences of internal variability—these snapshot periods nonetheless include a wide range of potentially relevant internal ocean-atmospheric oscillations.

We also note that although real-world greenhouse gas forcings are likely to be lower than assumed in the RCP8.5 scenario (), this is the only scenario for which high-frequency data are available as part of the CESM1-LENS dataset (). We further emphasize that although RCP8.5 is considered to be a high warming scenario, we explicitly intend to design a plausible “worst case scenario” storm and flood sequence in this analysis, and therefore, the use of a high-end emissions trajectory is appropriate.

We then rank all such 30-day cumulative precipitation events from each CESM1-LENS snapshot period, drawing from an effective sample size of 400-model years in each instance (10 years × 40 ensemble members). To ensure statistical independence of the dataset and that long-lasting events are not double counted, we require at least a 30-day separation between storm sequences. From among the top 3 ranked events in each period, we manually select a single 30-day storm sequence that exhibits large precipitation intensity peaks in both northern and southern California, as well as a pattern of 30-day cumulative precipitation that is spatially well distributed throughout both northern and southern portions of the state. This subjective aspect of the extreme event scenario selection process is critically important from the broader perspective of ARkStorm 2.0, which is designed to be a statewide exercise in which flood and emergency management capacity is severely tested. Therefore, we manually selected the respective ARkHist and ARkFuture events from among the top three ranked events such that each would bring a high level of impacts to the entire state rather than just a portion of the region. In so doing, we ultimately select the second ranked event for ARkHist (calendar date range: 2 September 2002 to 3 December 2002 in ensemble member #20) and the third ranked event for ARkFuture (calendar date range: 11 January 2072 to 11 February 2072 in ensemble member #2). Further analysis suggests that the selected ARkHist event has an approximate RI of ~85 years in the 1971–2020 era climate, and the ARkFuture event has an approximate RI of ~333 years in a 2051–2100 era high warming climate and is empirically unprecedented (i.e., a >400-year RI) in the 1971–2020 era climate (fig. S11).

LENS-WRF event-targeted downscaling approach

For each selected 30-day storm sequence, we use a high-resolution (3 km), nonhydrostatic regional weather model (WRF V4.3) embedded within initial and boundary conditions from CESM1 large ensemble (a framework known as “LENS-WRF”) to perform dynamical downscaling as originally developed by (). We use a full suite of three-dimensional atmospheric initial and boundary conditions from the high-frequency (6 hourly) temporal data available from the CESM1-LENS output files and conduct ~50-day long WRF simulations for each 30-day scenario event (allowing for ~1 week of model spin-up and ~1 week of event follow-up). Land surface initial and boundary conditions (including three-dimensional soil temperature, soil moisture, and snow depth) are drawn from the corresponding model member at monthly frequency (as this is the highest temporal resolution retained for three-dimensional land surface conditions in CESM1-LENS) such that they are spatiotemporally congruent with the atmospheric conditions.

In this analysis, we use a nonhydrostatic configuration of WRF-ARW (V4.3) including four nested domains with progressively finer spatial resolutions of 81, 27, 9, and 3 km (see fig. S12 for the detailed domain configuration). The outer three domains cover a large portion of the northeastern Pacific Ocean and the innermost 3-km domain also covers a broad oceanic region—as well as all of California and Nevada—to better represent near-coastal processes and sea-air interactions. WRF is configured using 44 vertical levels (with model top pressure at 50 hPa and vertical velocity damping turned on) and forced with time-varying SST (from CESM1-LENS). A higher density of vertical levels is prescribed near the surface to improve the representation of lower-level processes.

WRF physics parameterizations applied in these simulations include the Thompson graupel scheme (), the Kain-Fritsch (new Eta) cumulus scheme () (for 81-, 27-, and 9-km domains only; cumulus parameterizations are turned off for the innermost 3-km domain), the Dudhia shortwave radiation scheme (); the “rrtm” longwave radiation scheme (), the Yonsei University (YSU) boundary layer scheme (), the revised MM5 Monin-Obukhov surface layer scheme (), and the Noah-multiple parameterization (MP) land surface model (). The Noah-MP model includes a multilayer snowpack capable of liquid water storage and melt/refreeze cycles, direct representation of heat exchange due to phase changes, and a snow interception component allowing for canopy interception ().

Model validation and fitness for purpose

The overall performance of both CESM (as implemented in CESM1-LENS) and WRF have been previously assessed and validated in the context of both mean and extreme cool season precipitation in California. Swain et al. () found that the simulated distribution of CESM1-LENS cool-season precipitation was statistically indistinguishable from observations during the recent historical period in both northern and southern California. In addition, Huang et al. () found that high-resolution (3 km), nonhydrostatic WRF simulations nested within boundary and initial conditions from atmospheric reanalysis (i.e., pseudo-observations) were capable of simulating real-world extreme AR events (including extreme IVT) and associated extreme precipitation—including spatial patterns of orographic enhancement. However, we acknowledge that this validation does not obviate the potential for parametric and/or structural uncertainties that could lead to model biases that are difficult to quantify (as it is not possible to directly validate large ensemble climate model representation of specific extreme events). Nonetheless, the LENS-WRF configuration used in the present analysis is capable of generating physically realistic extreme storm events and is an appropriate tool for use in the context of “plausible worst case” scenario development.

Contextualization of CESM1-LENS relative to other large ensembles

We conduct additional analysis using daily precipitation data from several other large single-model ensembles [the 50-member CanESM2 (Canadian Earth System Model, Second Generation) at ~2.8° × 2.8° horizontal resolution, 20-member GFDL-CM3 (Geophysical Fluid Dynamics Laboratory Coupled Model, Version 3) at 2.0° × 2.5° horizontal resolution, and 30-member CSIRO-Mk3.6 (Commonwealth Scientific and Industrial Research Organisation Model, Version 3) at ~1.875° × 1.875° horizontal resolution] to aid in contextualization of the study’s primary focus on results driven by CESM1-LENS (40 members at 1° × 1° horizontal resolution). We note that CESM1-LENS has the highest horizontal resolution, by a wide margin, as well as the second largest number of ensemble members of these four large ensembles. To conduct as systematic an intercomparison as possible, we extract precipitation data for each of the top 4 ranked events in each ensemble and during each ARkHist and ARkFuture snapshot period. The results of this analysis are discussed in Results and can be visualized in figs. S1 and S2.

HUC region precipitation and runoff analysis

We select two “four-digit/subregional” HUC regions, as defined by the USGS, for more detailed analysis of regional precipitation and surface runoff during ARkHist and ARkFuture scenarios: HUC 1802 (Sacramento subregion, which includes the Sacramento River basin and Goose Lake watershed) and HUC 1804 (San Joaquin subregion, which includes the San Joaquin River basin; see fig. S10 for geographic outlines). We select these HUC regions, particularly, because they encompass most or all of the major SN western slope water storage and flood control reservoirs, as well as broad swaths of land in California’s Central Valley that are highly susceptible to large-scale flooding and are home to numerous flood control structures. We extract precipitation and runoff data from the WRF 3-km domain at 1 hour frequencies from geographic regions delineated by the respective HUC subregion shapefiles made available via the USGS (at https://apps.nationalmap.gov/downloader). We then plot empirical histograms of the upper tail of the precipitation (all values above 1 mm/hour) and runoff (all values above 0.5 mm/hour) distributions for each selected HUC region temporally aggregated at two different durations (1 and 24 hours) in both historical and future scenarios (Fig. 4).

Public availability of ARkStorm 2.0 atmospheric simulation data

Boundary and initial condition input files (derived from CESM1-LENS) and output files from the WRF simulations are archived on the DesignSafe web platform () via DOI: 10.17603/ds2-mzgn-cy51 ().