Atmospheric river activity during the late Holocene exceeds modern range of variability in California

Atmospheric rivers are associated with some of the largest ﬂ ood-producing precipitation events in western North America, particularly California. Insight into past extreme precipitation can be reconstructed from sedimentary archives on millennial timescales. Here we document atmospheric river activity near Leonard Lake, California, over 3,200 years, using a key metric of atmospheric river intensity, that is silicon/aluminum enriched layers that are highly correlated with modern records of integrated vapor transport. The late twentieth century had the highest median integrated vapor transport since the onset of the Medieval Climate Anomaly, with integrated vapor transport increasing during the Little Ice Age. The reconstruction suggests California has experienced pluvial episodes that exceeded any in the meteorologic instrumental era, with the largest episodes occurring two and three millennia ago. These results provide critical data to help avoid underestimation of potential risks and aid future planning scenarios

sediments of Leonard Lake; the model was then applied downcore.Our approach also incorporates time-uncertainty analysis, which is critical for quantitative interpretation of records in the paleosciences 28 .

Results and discussion
Extreme precipitation deposits Prominent light greenish-gray clay layers (10GY 7/1) 29 occur throughout the core, within a varying matrix of fine laminations and homogeneous zones (Fig. 2a, b).Finely laminated horizontal stratigraphy indicates deposition unimpacted by bioturbation.Downcore changes in grain size correspond to visually distinct stratigraphy such as clay layers and fine laminations (Fig. 2b, c).Fine laminations positively correspond with Fe/Al, which negatively correlates with homogenous zones (Fig. S1).In addition, high Si/Al values align with clay layers (Fig. 2d).Clay values above ~9% tend to produce a clay layer detectable to the naked eye and the correlation between Si/Al and percent clay is positive and significant (r = 0.30, p = 0.02, where autocorrelation was accounted for with the isospectral approach).Clay percentage data from the section of the core corresponding to the instrumental period encompass most of the range of percentage values present in the record from −550 CE to present; however, between −1050 to −550 CE, percent clay values exceed the highest values from the instrumental period.Measurements of Si/Al from the instrumental period encompass almost the entire range of variability over the full record; outside of the calibration period, only 2% of the values are higher than the Si/Al value in 1950 CE and none are lower than 1990 CE (Fig. 3a).Principal component analysis verifies trends in correlations, that is, grouping of siliciclastic elements (Si, Ti, Al, K); it also shows grouping of organic matter with Br on PC1 and the primary control of PC2 by Fe (Fig. S2, and additional sensitivity analysis in Fig. S3).

Correlation and proxy calibration
Of the proxies, Si/Al is most strongly correlated with integrated vapor transport (IVT, kg m −1 s −1 ) data (r median = 0.63, p median = 0.02), where all ensemble-member correlations are significant, and p values were corrected to account for false discoveries 30 (Fig. 3a, b; Fig. 4a).Strong, significant associations were also found between Si/Al and integrated water vapor (IVW; r median = 0.48, p median = 0.09) with significant ensemble-member correlations 96% of the time (Fig. 4b).Correlations between clay, IVT, and IWV are presented in the Supplementary Note 2 (Fig. S4a, b).Correlations between Si/Al and clay with Russian River (USGS station no.11461500, ~20 km from Leonard Lake) and Navarro River (USGS station no.11468000, ~30 km from Leonard Lake) streamflow data 31 were positive but not significant and are not considered further.Reconstructed IVT incorporates uncertainties in both the calibration regression statistics and the age model using an ensemble approach (Fig. 5a-e).The last occurrence of median IVT values as high as the modern period (~40,000 kg m −1 s −1 ) was before the onset of the Medieval warming in 950 CE, after which it appears to have declined for several centuries to the record's nadir of ~27,500 kg m −1 s −1 in the 1500 s (Fig. 5e).Median IVT rose with some fluctuation during the Little Ice Age (1300-1850 CE).The largest median peaks in IVT occur ca.−860 CE and −70 CE and were ~44,000 kg m −1 s −1 (Fig. 5e).

Extreme precipitation over the late Holocene in California
Stratigraphic data from Leonard Lake show remarkable repeating patterns throughout the core that we interpret as hydrologically forced changes in the lake's depositional processes.Prominent clay layers are interspersed with fine laminations and/or homogenous packages and, in the lower part of the core, occasional sandy layers (Fig. 2a, b).These lithological data suggest a process whereby extreme precipitation related to ARs transports fine siliciclastic rich material to the depocenter.The largest area of the watershed drains into the southeastern corner of the lake via a creek that traverses ~300 meters of meadow prior to entering open water.This creek overtops its channel and floods the meadow during high precipitation events and the coring site/depocenter lies ~500 m from its mouth.As a result, coarser Fig. 1 | Location and geomorphic features around Leonard Lake. a Leonard Lake is situated in the northwestern Coast Ranges (dark gray shading) 57 , California, USA 58 .b Streams (blue lines) 59 and the watershed (black line) are depicted around Leonard Lake (light blue) and its meadow (dark gray), with the world hillshade map layer 60 .
c Bathymetric map of Leonard Lake where increasing depth (m) is shown by increasingly dark blues, and the coring site is shown by a yellow circle in the lake's depocenter.Copyright © 2024 Esri and its licensors.All rights reserved.material is not carried to the depocenter, but rather the clay fraction flushes into the lake and slowly settles out, with visible layers indicating periods of extreme precipitation.We note a dominant pattern of homogenous sediment followed by a clay layer, overlain by a section of fine laminations.Fine laminations (≤mm scale) suggest suboxic to anoxic benthic conditions, allowing little degradation of organic matter and a lack of bioturbation.In addition, the correspondence between fine laminations and increased Fe/Al support the idea that sedimentary Fe was created in situ as a result of anaerobic decomposition (Fig. S1).This sequence suggests a process in which large storms raise lake levels, creating anoxic bottom water conditions that facilitate the preservation of fine laminations until (or unless) relatively dry conditions lowered lake levels.Dry conditions are likely recorded by homogenous packages that indicate oxygenated bottom waters and bioturbation.Modern lake level data from Leonard Lake support the underlying idea that large precipitation events appreciably raise lake levels while dry conditions lower lake levels (Supplementary Note 3, Fig. S5).The presence of sandy layers in the lower portion of the core may have formed concurrent with lower lake levels during nascent formation.The shallow lake levels would likely result in discharge much closer to the coring site because the mouth of the creek extended further into the basin.Consequently, coarser material could be more readily carried to the depocenter during high energy events.
Our evidence supports the premise that Leonard Lake consistently registered extreme precipitation over a 3200-year period.Faithful recording of extreme precipitation as clay deposition has allowed a quantitative reconstruction of AR-related metrics in the late Holocene, the first of its kind for the western United States.The reconstructed median IVT suggests the most recent analogue to modern extreme precipitation occurred over a millennia ago and that since 1300 CE there has been an upward trend.As a key metric of ARs, the IVT reconstruction indicates California has been subject to pluvial episodes that exceed any in the meteorological instrumental era, with the largest episodes occurring two and three millennia ago.
Based on inferences from proxy records, Northern Hemisphere midlatitude net precipitation has increased since 8 ka 32 and midlatitude Northern hemisphere models show storm activity has increased since 10 ka 33 .The early to middle Holocene was drier and a potential driver of these conditions was a weaker latitudinal temperature gradient.A weaker latitudinal gradient has again developed under anthropogenic warming; it has been argued that this could be consistent with future aridity and a reduction in mid-latitude water resources 32 .Our reconstruction appears to show that AR activity may not be linked to the northern hemisphere temperature gradient.In contrast, we find a slight decreasing trend in IVT values between −1300 and 1500 CE (Fig. 5e) while the latitudinal temperature gradient was strong, and an increasing trend in IVT values between 1500 CE and the present day (Fig. 5e) as the latitudinal temperature gradient weakened.The difference between our record's trends and others' 32 may simply reflect the fact that ARs are distinct features of the hydroclimate that are not interchangeable with net precipitation.Differences could also be related to sub-regional variability detected by our site and/or changes in temperature and insolation gradients, such as the role of cooler tropical Pacific sea surface temperatures and North American aridity 34 .Nevertheless, our work is consistent with recent model output comparing midlatitude landfalling ARs in western North America between the mid-Holocene and the preindustrial period 35 .This research shows a drying trend in annual AR precipitation between the two periods, which is supported by the decreasing IVT trend in our record before the preindustrial era, and is attributed to decreased atmospheric vapor content related to decreasing summer insolation 35 .Thermodynamic considerations, such as increased AR vapor content due to warmer air and sea surface temperature, are expected to predominantly drive increases in the number of ARs in the future 36 .Mechanism(s) notwithstanding, these data provide a critical perspective to improve the precision and accuracy of risk assessments for precipitation related threats, such as enhanced flash flooding or debris flow 13 .
Leonard Lake data align with regional climate reconstructions and suggest the reliability of Si/Al to record both periods of increased precipitation as well as drought (Fig. S6a-c).For example, we find increasing IVT during the Little Ice Age consistent with cooler and wetter climatic expectations.Heavy precipitation, as reconstructed from blue oak tree rings in northern California, also suggest an increase in the frequency or intensity of ARs in northern California since the 1700s 37 (Fig. S6a).Our reconstruction also indicates decreasing IVT during the Medieval Climate Anomaly, a period characterized by centennial-scale aridity in the North American West 37 .Our reconstruction shows several eras of low IVT occurring in 600-725 and 1250-1350 CE, as well as the 1500 s-the latter of which is concurrent with the widespread 16th century mega-drought 38,39 .A reconstruction of Sacramento River flows found the driest 20-year period was in the early 1300s 40 , aligning with the low IVT reconstructed in this study between 1250-1350 CE.The Sacramento River flow record 40 also shows some of the driest individual years occurred in the early 1500 s (Fig. S6b).The consensus across multiple records lends credence to our reconstruction.Although our record cannot reconstruct specific storm  events, the alignment with these broader climate trends highlights the importance of extreme precipitation in driving pluvial/drought cycles in California going back millennia.
Although the uncertainty analysis in this work considered a range of Si/ Al-AR linkages, we cannot account for the possibility that the relationship between Si/Al and extreme precipitation changed through time (non-stationarity), a problem besetting all paleo-reconstructions.Other variables that may affect siliciclastic transport and preservation within the sedimentary archive include the changes in lake level mentioned above and vegetation changes that alter run-off to the site.However, we suggest that the processes influencing Si/Al and clay content at Leonard Lake have been sufficiently stable to support a consistent proxy signal for millennia.Leonard Lake is situated in the ancestral territory of the Pomo Tribe.Tribal members call this lake Ka-su'-su, meaning waters malevolent or waters bedeviled, likely due to the extremely steep shore and surprising depth; Indigenous people are thought to have avoided using the lake and surrounding area 41 .Hence, there is evidence the lake has long been deep and the area relatively undisturbed.In the recent past, the old-growth redwood forests in Leonard Lake's watershed have been protected by the landowners, unlike the surrounding area which experienced heavy logging 41 .The relatively constant sedimentation rate and lack of obvious major sedimentological changes also supports the idea of long-term stability.
The recognition that a paleo perspective is necessary to accurately understand the risks associated with extreme events 13 continues to motivate research concerning the long-term hydrologic history of the western United States.Evidence from this study suggests that paleo records of AR activity are detectable given the right long-term depositional environment, such as is found at Leonard Lake.Our reconstruction corresponds with recognized paleoclimatic periods including the Medieval Climate Anomaly and Little Ice Age, corroborates regional tree ring reconstructions, and is consistent with a mega-drought event.Importantly, these findings show that extreme precipitation events have been a mechanism driving both pluvial and drought conditions in California for at least 3200 years.This study contextualizes present day storms by demonstrating the range of AR-related extreme precipitation on millennial time scales in the hydrologically dynamic area of California.

Methods
In April 2014, two surface cores and three sets of overlapping long cores were recovered from the depocenter of Leonard Lake.Bolivia and Livingston corers were used to obtain overlapping cores with a total length of 446 cm.Core liners were used to preserve stratigraphic integrity and to aid core transport and sampling.All cores were split, described, imaged, and sub-sampled in the Quaternary Paleoenvironmental Research Laboratory at the U.S. Geological Survey in Menlo Park, California.Lake depth measurements were collected with a Norcross H22PX® and geolocated with a Trimble GeoXH® to create a bathymetric map (Fig. 1c).

Geochronology
A total of 17 samples of discrete organic materials were obtained for AMS 14 C dating and were processed at the National Ocean Sciences Accelerator Mass Spectrometry center at Woods Hole Oceanographic Institution, the US Geological Survey Radiocarbon Laboratory, and Beta Analytic (Table S1).We note the date at 236 cm represents an age reversal and was treated as an outlier by the model.Given that the two proximal dates at 237.5 cm and 244.5 cm overlap at an uncertainty of two sigma, we deem the exclusion of the date at 236 cm as justifiable.Short-lived radionuclides ( 210 Pb, 137 Cs, and 226 Ra) from the upper 60 cm were measured using gamma spectrometry 42 (Table S2).Seventeen samples were measured until unsupported 210 Pb was no longer detected and a definitive peak 137 Cs was found.Counting times of 48 to 72 h were generally required to reach a statistical error of less than 10% for excess 210 Pb in the deepest samples.Excess 210 Pb was calculated as the difference between total 210 Pb and 226 Ra activities.Downcore 137 Cs values show a peak at 28 cm (Table S2), indicating fallout from nuclear weapons testing in 1963 CE. Background 210 Pb was likely reached around 50 cm, given the low and static dpm/g of the deepest 210 Pb measurements (Table S2).These data were used to construct the agedepth model (Fig. S7a, b) using the Bayesian software Plum (version 0.3.0 using default priors) 43 , in RStudio 44 .In Plum, Monte Carlo age-depth fitting uses calibrated-age probability distributions, enabling uncertainty calculations at every level.
To understand the effect of age model uncertainty on our results, the geoChronR package in R was used to propagate error and visualize results 28 .Radiometric data and proxy data were reformatted into a Linked PaleoData (LiPD) dataset.GeoChronR requires users to run one of four age-modeling techniques, of which the algorithm in Bacon most closely aligns with our approach 45,46 (Plum is not available in geoChronR).The underlying process is as follows: within geoChronR, Bacon is run normally and returns a subset of the Markov Chain Monte Carlo accumulation rate ensemble members with high a posteriori probabilities; then, an age ensemble-a composite of 1000 highly probable age models-is extracted for analyses.To incorporate 210 Pb error in our analyses, we ran Bacon with seventeen modeled 210 Pb ages and their associated uncertainties from Plum (Fig. S8, Table S3) as well as 14 C data.Both the Plum and Bacon models produced highly similar age estimates (average of 3.3% difference; Fig. S9) for the full age models, justifying this approach.The mean deposition time over the entire core is 7.4 year cm −1 with a mean of 1.98 year cm −1 during the calibration period (1948-2015 CE).Average model uncertainty is ±74.5 year over the 3200-year-long core, while average uncertainty through the calibration period is ±3.5 year.

Stratigraphic, grain size, and geochemical analyses
The sediment was described, at 0.5-cm resolution, as one of four types: clay layers, fine laminations, homogenous deposition, and sand layers.Grain size measurements were carried out contiguously at 0.5-cm resolution.Samples were processed first with 50 mL of bleach (5.25% sodium hypochlorite) for 24 h to remove organic material, water washed, and then treated with 10% HCl to remove carbonates.Analysis was conducted on a Beckman Coulter LS13 320 Laser Diffraction Particle Size Analyzer in the Unsaturated Flow Processes laboratory at the U.S. Geological Survey.Water content at 50 °C, loss on ignition (LOI) at 550 °C (organic matter), and at 960 °C (carbonate) were determined at 1-cm contiguous intervals 47 ; the residual is quantified as the non-carbonate inorganic fraction.A Thermo Scientific Niton XL3t GOLDD portable XRF (p-XRF) was used to collect 1-cm resolution elemental data on split cores with standard reference materials: 165a 48 , 2709a 49 , and JMS-1 marine sediment 50 .

Analyses and instrumental data
Principal component analysis (PCA) was undertaken on elemental and LOI datasets using compositional data analysis techniques [51][52][53] (see Supplementary Methods and Supplementary Note 1).Pearson's correlation coefficients between proxy and instrumental data were calculated for each age ensemble 28,54 .The false discovery rate was controlled 30 .We calibrated the instrumental record of ARs with values of Si/Al and percent clay using age assignments derived from the age-depth model for each sample depth and regressing these against the dataset of AR metrics 25 , following established protocol 26 .
A key feature of ARs is that they move large amounts of water vapor.Hence, metrics that incorporate moisture and wind best reflect this water vapor flux 25 .ARs have been defined by vertically integrated horizontal vapor transport (IVT, kg m −1 s −1 ) and integrated water vapor (IWV, kg m −2 ) exceeding defined thresholds 5 .As such, we calculated an annual sum for both IVT and IWV for ARs that were specified at 40°N in Gershunov's 25 database 55 .We then calibrated our proxy data to these sums.For each iteration in the analysis, an age ensemble-member was selected and used to bin the Si/Al (or clay data) onto a 5-year interval where data were averaged.The instrumental AR data were also binned at the same timescale.Streamflow data (ft 3 sec −1 ) from stream gages in two local rivers, the Russian River (USGS station no.11461500) and Navarro River (USGS station no.11468000) 31 , were summed annually between 1950-2014 and considered along with the instrumental AR data.Within geoChronR, a linear regression model was derived, and then proxy data were scaled to guard against slope attenuation ("regression dilution") 56 .The model was used to reconstruct integrated vapor transport (IVT sum kg m −1 s −1 ) from proxy data back in time, with age uncertainty propagated through the regression.

Fig. 2 |
Fig. 2 | Lithological and geochemical results from Leonard Lake.Core analyses from Leonard Lake are plotted on the same depth scale (left): a Composite line scan imagery; b lithological descriptions based on visual analysis where clay layers are blue, fine laminations are black, homogenous areas are yellow, and sand layers are red and; c percentage data for clay (blue), total sand (gray), silt (orange), and noncarbonate inorganics (yellow); and (d) Si/Al values (black).

Fig. 3 |
Fig. 3 | The proxy data (Si/Al) and instrumental data (IVT) during the calibration period.a The graph shows variation in Si/Al with age uncertainty bands; the median estimate is shown as a solid black line and the 50 and 95% highest probability density regions are shown in dark and light gray.b The graph shows the calibration-in-time predictand where the y-axis is annually summed IVT (kg m −1 s −1 ) and binned at a 5-year interval.

Fig. 4 |
Fig. 4 | Correlation distributions of Si/Al and IVT and Si/Al and IWV and their significance.a The distribution of Si/Al and IVT and (b) the distribution of Si/Al and IWV using the "effective-n" approach 28,53 with the false discovery rate (FDR) controlled.The ensemble count is on the y-axis and the r value is on the x-axis.In (a), the median r value (0.63) is shown by the vertical red line and the 0.05 and 0.95 percentile r values are shown by dashed red lines.In (b), the median r value (0.48) is shown by the vertical red line and the 0.05 and 0.95 percentile r values are shown by dashed red lines.

Fig. 5 |
Fig. 5 | Calibrated reconstruction of paleo-IVT using ensemble regression.a The graph shows the calibration-in-time predictand.b The scatterplot shows the relation between Si/Al and IVT with points from 100 age ensemble members plotted with transparency and ensemble models plotted as partially transparent red lines.The right panel shows the distribution of the regression model (c) slope and (d) shows the intercepts, where vertical red lines indicate the 50th (solid), 2.5th, and 97.5th (dashed) percentiles.e The graph shows reconstructed IVT using age ensembles and regression model ensembles.The median estimate is shown as a solid black line and the 50 and 95% highest probability density regions are shown in dark and light gray.