Holocene hydroclimatic variability recorded in sediments from Maddox Lake (northern California Coast Range)

Abstract

Perspectives on past climate using lake sediments are critical for assessing modern and future climate change. These perspectives are especially important for water-stressed regions such as the western United States. One such region is northwestern California (CA), where Holocene-length hydroclimatic records are scarce. Here, we present a 9000-year, relative lake level record from Maddox Lake (CA) using a multi-indicator approach. The Early Holocene is characterized by variably low lake levels with a brief excursion to wetter climates/relative highstand ca. 8.4–8.06 cal ka BP, possibly related to the 8.2 ka cold event and changing Atlantic Meridional Overturning Circulation (AMOC). From 5.2–0.55 cal ka BP, Maddox Lake experienced a long-term regression, tracking changes in summer-winter insolation, tropical and northeast Pacific SSTs, and the southward migration of the ITCZ. This gradual regression culminated in a pronounced relative lowstand during the Medieval Climatic Anomaly (MCA). A marked relative highstand followed the MCA, correlative to the Little Ice Age. The latter reflects a far-field response to North Atlantic volcanism, solar variability, and possibly changes in AMOC and Arctic sea ice extent. Our results further confirm the hydroclimatic sensitivity of northwest California to various forcings including those emanating from the North Atlantic.

INTRODUCTION

California faces a perennial freshwater availability crisis (MacDonald et al., 2008; Diffenbaugh et al., 2015; Hatchett and McEvoy, 2018). Key to the crisis is the feast-or-famine, winter-dominated climatology (Wang et al., 2017; Swain et al., 2018). Climate models suggest that future CA climate change will be characterized by an increase in the magnitude of droughts and pluvials, as well as enhanced “whiplash” drought-to-pluvial winter climate response (Das et al., 2013; Swain et al., 2018; Ullrich et al., 2018). Critical to understanding this future projection is a clear view of natural climate variability and its climatic forcings.

To accomplish this, paleoperspectives are required to provide a baseline understanding of natural climate variability and magnitude as well as their drivers (i.e., climatic forcings) that modulate climate at annual to millennial timescales (Schmidt et al., 2004; Shuman and Marsicek, 2016; Bakker et al., 2017; Cook et al., 2018; Lachniet et al., 2020; Tierney et al., 2020). These paleoperspectives are especially important for key geographical regions in northwestern California (CA) where changes in water availability affect water/reservoir storage/management and water distribution to more arid CA regions. One such region is the northern California Coast Range and Klamath Mountains. Two of California's largest reservoirs (Shasta and Trinity) lie within this key region of water replenishment. Yet, our understanding of Holocene hydroclimatic variability using analyses similar to those presented in this study remains underdeveloped for this important region. To improve this understanding, we are investigating regional lakes as natural barometers of winter hydroclimatic variability.

Here, we present a new 9000-year record of relative lake level change (Maddox Lake, CA) in the northern California Coast Range (Fig. 1). We use a combination of physical, chemical, and biological indicators to identify the dominant sediment-climate signal. Our interpretations include variably low lake levels during the Early Holocene interrupted by a relative highstand ca. 8.4–8.06 cal ka BP, possibly related to the 8.2 ka cold event, a Middle Holocene regression into a prominent Medieval Climatic Anomaly (MCA) relative lowstand, and a short-lived Little Ice Age (LIA) relative highstand. Overall, the past 9000 years were dominated by lake level regression, reflecting the evolution of winter–summer insolation, migration of the intertropical convergence zone (ITCZ), long-term changes in Pacific Ocean-atmosphere forcing, and their combined modulation of the magnitude and frequency of winter storm tracks across CA. There are, however, notable excursions from this low frequency signal including the ca. 8.4–8.06 cal ka highstand (possibly related to the 8.2 ka cold event) and the LIA highstand. The latter two events represent departures from the dominant low frequency Holocene forcings and suggest far-field responses to climate drivers such as Atlantic Meridional Overturning Circulation, volcanism, and solar variability.

Figure 1. Study site location map and regional perspective. Regional sites mentioned in the text: 1. ODP Site 1019 (Barron et al., 2003), 2. TN062-O550 (Barron et al., 2018), 3. Oregon Caves National Monument (Ersek et al., 2012), 4. Upper Squaw Lake (OR) (Colombaroli and Gavin, 2010), 5. Sanger Lake (Briles et al., 2008, 2011; Briles, 2017), 6. Bolan Lake (OR) (Briles et al., 2005, 2008; Whitlock et al., 2008), 7. Twin Lakes (Wanket, 2002) and Fish Lake (Crawford et al., 2015), 8. Lake Ogaromtoc (Crawford et al., 2015), 9. Campbell Lake (Briles et al., 2011; Briles, 2017) and Taylor Lake (Briles et al., 2011; Briles, 2017), 10. Crater Lake (CA) (Mohr et al., 2000), Bluff Lake (Mohr et al., 2000), Cedar Lake (Briles et al., 2011; Briles, 2017), and Mumbo Lake (Daniels et al., 2005), 11. Flycatcher Basin (R.S. Anderson et al., 2008), 12. White Moon Cave (Oster et al., 2017).

BACKGROUND

Regional climatology

The northern California Coast Range is characterized by a Mediterranean climate with wet, cold winters and warm, dry summers. This winter-dominated climatology is controlled by the position of the winter polar front in relation to the Pacific subtropical high (Cayan and Peterson, 1989; Cayan et al., 1998; Dettinger et al., 1998; Wise, 2010). In general, a weak winter season Pacific subtropical high corresponds to wetter winters in California as storms track more frequently across the state. Notably, atmospheric rivers account for most of this winter season moisture transport into CA, generating alpine snowfall, lower-elevation precipitation, and occasional flooding (Ralph et al., 2006; Dettinger et al., 2011).

Superimposed on this unimodal winter climatology is the western United States precipitation dipole (Cayan et al., 1998; Dettinger et al., 1998; Fye et al., 2004; Wise, 2010, 2016). The dipole is characterized as a transition zone separating a N-S antiphased precipitation regime. Based on a 500-yr tree ring dipole reconstruction, average dipole latitude is 40°N (i.e., the location of Maddox Lake); although, the dipole is highly variable over time, ranging from 35–44°N over the 500 yr period of study (Wise, 2010, 2016). Of course, this boundary is dynamic and occasionally absent altogether (Cayan et al., 1998; Dettinger et al., 1998; Fye et al., 2004; Wise, 2010, 2016). As the latitude of the dipole changes from year to year, the mean land-falling latitude of winter storm tracks also changes. Furthermore, it is well documented over the historical period that relatively small interannual changes in the frequency and/or magnitude of winter storms can generate large changes in California lake and reservoir hydrology (Benson et al., 2002; Kirby et al., 2004, 2006; Hanson et al., 2006; Adams et al., 2015).

Wise (2010) also highlighted the dipole's relationship from 1926–2007 AD to larger-scale ocean-atmosphere dynamics such as El Niño-Southern Oscillation (ENSO/SOI), the Pacific Decadal Oscillation (PDO), and the Atlantic Multidecadal Oscillation (AMO). In terms of the Pacific, over CA, the dipole is most prominent during −PDO/+SOI (dry south/wet north) and +PDO/−SOI (wet south/dry north) phases and weakest during −PDO/−SOI and +SOI/+PDO phases (Wise, 2010). Thus, the dipole is particularly sensitive to SOI variability and conditions in the tropical Pacific as well as the extra-tropical Pacific (Cayan et al., 1998, 1999; Dettinger et al., 1998; Wise, 2010; Peng et al., 2013). In fact, model analyses attribute up to a third of CA winter precipitation variability between 1895–2014 AD to Pacific SST forcing (Seager and Hoerling, 2014). Of course, the role the tropics play in CA's hydroclimate varies over time in response to internal variability, changing climatic boundary conditions, and/or other far-field ocean-atmosphere dynamics (Seager et al., 2005, 2014; Cook et al., 2011, 2013; Chiang and Friedman, 2012; Peng et al., 2013; Coats et al., 2016; Kam and Sheffield, 2016; Wise, 2016). In the Atlantic, the AMO is an important player in modulating North Atlantic climate. As Wise (2010) and Zhang and Delworth (2007) demonstrated, the AMO affects Pacific climate as well. Wise (2010) showed that the −AMO, when concurrent with −SOI or +SOI creates a prominent dipole structure, very similar to the relationship between the PDO and SOI.

Study site

Maddox Lake is landslide formed and sits at 857 m above sea level in Trinity County, CA (Fig. 1). The surface area of Maddox Lake is ~0.04 km² with a drainage basin of ~1.54 km². The main inlets enter from the south and northwest, although no inlets were active as of July 2019. There is an outlet at the east end of the lake (~1 m above modern lake water surface elevation); however, the outlet was not active at the time of core acquisition. Maddox Lake is uniformly shallow with emergent macrophytes (e.g., cattail [Typha latifolia]) dominating the extensive littoral zone, with bordering sedges (Cyperaceae) and rushes (Juncaceae). The surrounding understory was open with Adler (Alnus rhombifolia), California ash (Fraxinus dipetala), Douglas fir (Pseudotsuga menziesii), Pacific Madrone (Arbutus menziesii), California hazelnuts (Corylus cornuta), and Coast Live Oak (Quercus agrifolia). As a remote, isolated site, there are no historical limnological or surface-elevation data from Maddox Lake. Using Google Earth historical imagery, there are only 14 interpretable (i.e., in focus) images from 1985–2021, and only two images prior to 2004. With the exception of seasonal changes in vegetation in and around the lake, changes in lake level are not obvious with the images available. Water chemistry from the modern lake is characterized by pH (6.83 ± 0.05, n = 3), conductivity (214 μs ± 11.53, n = 3), temperature (21.27°C ± 0.38, n = 3), total dissolved solids (152 ppm ± 6.56, n = 3), and salinity (0.11‰ ± 0.01, n = 3).

Modern climate data were taken from Forest Glen, CA (22 km south of Maddox Lake and 693 m asl) for the period 1930–1985 AD. Precipitation averages 154 cm/yr, with >80% total precipitation between November and March; snowfall averages 67 cm/yr. Winter season (DJF) temperatures average 3.2°C; whereas summer (JJA) temperatures average 19°C. These data show that winter hydrology (liquid and solid) and subsequent snowmelt runoff control the lake's inputs, and summer evaporation likely is the major control on water loss. The extent that groundwater and direct water loss via the lake's outlet modulate the lake's hydrology is unknown.

Existing regional Holocene research

From the northern California Coast Range and surrounding region, there exist a variety of Holocene paleorecords including lacustrine, speleothem, marine, and tree-rings (Fig. 1) (Adam and West, 1983; West, 1993; Heusser, 1998; Mohr et al., 2000; Meko et al., 2001; Barron et al., 2003; Daniels et al., 2005; Vacco et al., 2005; Briles et al., 2008, 2011; Whitlock et al., 2008; Colombaroli and Gavin, 2010; Ersek et al., 2012; Malevich et al., 2013; Crawford et al., 2015; Briles, 2017; Anderson et al., 2020). The lacustrine records are predominantly fire (i.e., charcoal) and vegetation (i.e., pollen) reconstructions. In general, they show a pronounced warming and drying from the late glacial into the Early Holocene. Within the Holocene, the story is more complex due to site-to-site differences in topography, slope aspect, bedrock (soil) type, and distance (i.e., moisture gradient) from the Pacific Ocean (Briles, 2017). Despite some of these complex responses, there emerges a general picture of Holocene fire and vegetation dynamics.

Briles (2017) compiled charcoal and pollen data from seven lakes in the Klamath Mountains (northwest CA) to characterize Holocene vegetation and fire activity (Fig. 1). Fire activity was a persistent feature of the Holocene with some regional differences. Low-frequency changes, however, in fire severity seem to track summer insolation. Pollen data from the seven study lakes suggest that the Early Holocene was warmer and drier than present, the Middle Holocene was cooler and wetter than the Early Holocene, and the Late Holocene was drier than the Middle Holocene and drier than present (Briles, 2017). A closer look at the Late Holocene from two of these seven study lakes (Fish Lake and Lake Ogaromtoc) by Crawford et al. (2015) suggests a warmer and drier MCA and a cooler and wetter LIA for the region (Fig. 1). However, anthropogenic (i.e., Indigenous Peoples) influences on the vegetative structure of these lakes’ drainage basins warrant caution when interpreting paleoecological records in areas affected by human activity. Pollen and charcoal data from nearby Upper Squaw Lake (southern Oregon) suggest similar regional responses to the Late Holocene MCA and LIA (Fig. 1) (Colombaroli and Gavin, 2010).

The highest resolution (ca. 3 years/sample) Holocene hydroclimatic record from the region is found at Oregon Caves National Monument (Fig. 1) (Ersek et al., 2012). Using stable isotopes δ¹³C_{(speleothem calcite)} and δ¹⁸O_{(speleothem calcite)}, Ersek et al. (2012) inferred winter precipitation amount and/or surface biomass development and atmospheric temperature, respectively. Low-frequency Holocene δ¹⁸O_{(speleothem calcite)} changes likely are associated with winter insolation forcing and suggest a small (~1.0°C) increase in winter temperature from the Early to the Late Holocene; however, there is no similar low-frequency trend observed in the δ¹³C_{(speleothem calcite)}. Both isotopes also reveal multi-decadal to millennial scale variability, suggesting pronounced changes in both precipitation and temperature throughout the Holocene. Ersek et al. (2012) attributed this higher frequency change to a combination of solar forcing and Pacific Ocean-atmosphere dynamics.

From the adjacent northeast Pacific, Barron et al. (2003, 2018) combined micropaleontological data, reconstructed sea surface temperatures (SST), and pollen from two core sites (ODP site 1019 and TN062-O550) to characterize and explain changes in the marine environment and coastal vegetation (Fig. 1). In general, at ODP site 1019, SSTs were higher than modern in the Early Holocene, lower than modern in the Middle Holocene (8.2–3.2 cal ka BP), and higher than modern in the Late Holocene (Barron et al., 2003). Based on the pollen data, the coastal climate was characterized by warmer and drier conditions between 9–5.2 cal ka BP with generally warming winters, cool summers, and mild winters between 3.5 cal ka BP and modern. The north coast temperate rain forest was established ca. 5.2 cal ka BP, based on alder and coastal redwood pollen (Barron et al., 2003).

Located 96 km south of ODP site 1019 and 103 km northwest of Maddox Lake, marine core TN062-O550 indicates a long-term warming of winters from 7.3 cal ka BP to modern, punctuated by an abrupt shift in terrestrial (warmer winters) and marine (warmer SSTs) conditions ca. 2.8–2.6 cal ka BP (Fig. 1) (Barron et al., 2018). Between 4.2–2.8 cal ka BP, coastal winter precipitation potentially increased as well as fluvial discharge; thereafter, pronounced warm-cold cycles characterized the past 2.8 cal ka BP. Changes in insolation and Pacific Ocean-atmosphere dynamics linked to the Pacific Decadal Oscillation and El Niño-Southern Oscillation are considered the dominant drivers of Holocene variability at these marine sites (Barron et al., 2003, 2018; Barron and Anderson, 2011).

METHODS

Due to limited site access, a single multidrive sediment core was collected using a Russian coring system 15 m from the modern edge of Maddox Lake (40°33.210′N, 123°25.190′W) in 0.20 m water (MLRC18-1 [48–261 cm total core length over 3 drives]), and ~130 m north of the modern inlet (not flowing in 2018). This location was selected to capture lake level changes as imparted on the near shore sedimentology (Harrison and Digerfeldt, 1993; Dearing, 1997; Shuman and Serravezza, 2017). The core was opened on site, photographed, described, and wrapped in plastic wrap and aluminum foil for transport back to California State University Fullerton.

Age control for this study was determined using radiocarbon dating (accelerator mass spectrometry [AMS] ¹⁴C dating) of 17 discrete microscopic and macroscopic (>125 μm) organic materials (e.g., charcoal, pine needles, seeds). One sample (24.5 cm) returned a modern age and was not used in the age model (see results). Samples were measured on the insoluble fraction for ¹⁴C at the University of California, Irvine Keck Carbon Cycle AMS Facility (Table 1). All samples were treated with an acid-base-acid protocol (1N HCl and 1N NaOH, 75°C) prior to combustion. Radiocarbon concentrations are given as fractions of the modern standard, D¹⁴C, and conventional radiocarbon age, following the conventions of Stuiver and Polach (1977, p. 355). Sample preparation backgrounds have been subtracted, based on measurements of ¹⁴C-free wood. All results have been corrected for isotopic fractionation according to the conventions of Stuiver and Polach (1977), with δ¹³C values measured on prepared graphite using the AMS spectrometer (these can differ from δ¹³C of the original material and are not shown). Samples labeled “Modern” contain excess ¹⁴C, probably from the mid-twentieth century atmospheric thermonuclear weapons testing (Table 1).

Table 1. Radiocarbon data for cores MLRC18-1.

Magnetic susceptibility (MS) (× 10^-7 m³/kg) was measured at 1-cm contiguous intervals using a Bartington MS2 magnetic susceptibility meter. The same samples were used to determine percent water content and percent dry bulk density. Once dried, the samples were combusted at 550°C and 950°C for two hours to calculate the percent total organic matter (%TOM) and percent total carbonate (%TC), respectively (Dean, 1974). Mass accumulation rates were calculated by dividing the sedimentation rate per depth by the dry bulk density per depth (g/yr/cm²).

X-ray fluorescence (XRF) was measured at 1-cm intervals on the core surface with a plastic wrap surface barrier. We used a portable Olympus Vanta C-series pXRF with a 2 beam 40kV energy source in Geochem mode for 60 seconds to acquire elemental data. Measured data were calibrated using the Olympus calibration standard and exported in parts per million (PPM) using the equipment's internal software. Mo has a limit of detection of <5 ppm and Ti <25 ppm. A handheld XRF method for determining elemental concentrations in lacustrine sediments is a commonly used, rapid method (e.g., Niederman et al., 2021; Wright et al., 2023). However, the data should be considered semi-quantitative due to inherent sediment heterogeneity and its variable clastic, organic, and chemical components (Weltje and Tjallingii, 2008). As a result, we interpret the Ti and Mo data in the context of our other sediment measurements and their various interrelationships (e.g., MS, %TOM, grain size, etc.).

Oogonia capsules per 1 g dry sediment >125 μm were counted at 2-cm intervals (e.g., 0–1 cm = 0.5 cm, 2–3 cm = 2.5 cm, etc.). Oogonia are the reproductive capsules of the aquatic macrophyte Chara (sp.) (Groves and Bullock-Webster, 1924; Burne et al., 1980; Vance et al., 1992; Mullins, 1998; Apolinarska and Hammarlund, 2009; Détriché et al., 2009; Kemp et al., 2012). Only oogonia capsules with >50% remaining “body” were counted. If there were no oogonia >50% body preservation in the sample, but they were present in <50% parts only, a value of 1 oogonia was assigned for the depth. Notably, none of the oogonia contained their carbonate exterior. In each case, only the black carbon interior was preserved, suggesting that post-depositional dissolution of carbonate occurred.

Grain size was measured at 1-cm contiguous intervals following standard pretreatment protocols: 30–50 mL of 30% H₂O₂, 10 mL of 1N HCl, and 10 mL of 1N NaOH (Leidelmeijer et al., 2021). Grain size was determined using a Malvern Mastersizer 2000 grain size analyzer attached to a Hydro 2000 G dispersion unit. A 10-second sonication proceeded each analysis in the Hydro 2000 G dispersion unit prior to analysis. A silica carbide polishing powder standard was run twice at the beginning of each day, once every 10 samples, and once at the end of every day to evaluate the equipment's analytical stability over time (n = 4127, average = 13.11 μm, standard deviation = 0.10 μm). All data are reported as volume percent and divided into 10 grain-size intervals according to the Wentworth scale (Wentworth, 1922) as well as d0.5 (0.5 = mean).

We utilized Principal Component Analysis (PCA) based on normalized data and Euclidean distance to explore the relationships among magnetic susceptibility, percent total organic matter, percent sand, percent clay, titanium, and molybdenum concentrations (n = 211). Percent silt was not included because it was the dominant grain size fraction and overwhelmed PCA output. Oogonia also were not included because they were not sampled at 1-cm contiguous intervals. To help confirm the statistical significance of the relationships between samples, we calculated a similarity profile (SIMPROF) permutation test. SIMPROF analysis enabled us to test for structure in multivariate data and returned a P-value to determine whether the multivariate structure manifest in a group of samples is more or less similar to each other than would be expected if the data were random and lacked structure (a more detailed review of this technique can be found in Clarke et al., 2008, and Somerfield and Clarke, 2013). All statistical analyses were conducting using PRIMER V7 (PRIMER-E, Plymouth, UK) statistical program.

A PCA scatterplot was generated using the first and second principal component eigenvectors. All standardized data were plotted on the PCA scatterplot with symbols that correspond to the groups determined by cluster analysis. The position of each data point on the scatterplot was influenced by how strongly each principal component influenced each individual data point, respectively. The lines plotted in the scatter plot show the direction wherein each variable, or indicator, is increasing. The scatterplot analysis reveals which indicator accounts for most of the variance in the data for each cluster group.

RESULTS

The bottom 12 cm (261–249 cm) is a dark brown sediment with faint layering (Fig. 2). From 249–221 cm, the sediment is again a mottled brown to grayish green color, with a distinct dark brown unit from 240–236 cm. From 221–202 cm, the sediment is dark brown with minor grayish-green mottling; the change at 202 cm is abrupt. The sediment is very gray between 202–195 cm. Between 195–188 cm, the sediment again is distinctly medium brown with less grayish-green material. At 188 cm, the sediment becomes more grayish green with less brown mottling to 146 cm. Another abrupt change occurs at 146 cm, where the sediment becomes a mottled brown to grayish green color down to 112 cm. Within the latter unit is a clean light green interval between 135–133 cm. A change to variably dark brown organic-rich sediment occurs at 112 cm followed by an abrupt transition at 75 cm. The latter unit is characterized by faint layering of variegated brown mud. From 75 to 64 cm, the sediment is a light gray brown. From 64–48 cm, the sediment is massive, dark brown with some visible organics. The upper 48 cm of the core stratigraphy is a massive brown peat with abundant, visible roots. This upper section of the core is not included in the study because dense roots prohibited cm-scale sampling for a full suite of analyses. However, surface samples were measured for MS, %TOM, Ti, and Mo to characterize the modern lake sedimentology (see Discussion).

Figure 2. Core MLRC18-1 sediment data versus depth with location of calibrated ¹⁴C ages (far left side; values given in cal yr BP). From left to right: (A) Magnetic susceptibility (× 10^-7 m³/kg), (B) number of oogonia capsules per 1 g dry sediment >125 μm (C) percent total organic matter, (D) Mo concentration (ppm), (E) Ti concentration (ppm), (F) percent clay, (G) percent silt, (H) percent sand. Visual stratigraphic column is shown at the far right.

An age model was developed using 16 of the 17 radiocarbon dates (Fig. 3; Table 1). A “modern” ¹⁴C age at 24.5 cm was not included in the construction of the age model. The surface age at 0 cm was assumed 2018 AD (−68 calendar years before present [present = CE 1950]) for Bacon input. These dates were entered into the Bacon (v.2.3, IntCal13) age-modeling software, which incorporates underlying assumptions (e.g., sediment accumulation rates) into a Bayesian statistical model (Blaauw and Christen, 2011). The age model produced a median surface age of −65 yr BP and basal age of 8930 yr BP. Sedimentation rates average 0.08 cm/yr (± 0.06) [or 34.7 yr/cm (± 41.8)] with a maximum of 0.23 cm/yr [or 5.8 yrs./cm] and a minimum of 0.01 cm/yr [or 170.4 yr/cm]; Fig. 3).

Figure 3. (A) Maddox Lake core MLRC18-1 age–depth plot using the Bacon (v.2.2, IntCal13) age-modeling software (Blaauw and Christen, 2011). Blue features are the calibrated ¹⁴C dates; gray stippled lines show 95% confidence intervals. (B) X-axes show time in cal yr BP vs. y-axis, which shows sediment sample age resolution in yr/cm.

Magnetic susceptibility values range from −0.3 to 1.3 x 10^-7 m³/kg with an average of 0.4 ± 0.4 × 10^-7 m³/kg (Fig. 2). There are intervals of higher-than-average MS between 250–240 cm, 213–119 cm, and 73–65 cm. Oogonia range from 0–2308 counts/g dry sediment with an average of 213 ± 371 counts/g dry sediment; the results are plotted on a log scale (Fig. 2). In general, oogonia are present between 260–204 cm and 154–76 cm; however, there is significant variability throughout, with very few oogonia between 204–154 cm and 76–48 cm. Percent total organic matter ranges from 8.2–80.7 % with an average of 29.2 ± 18.7% (Fig. 2). The %TOM is uniformly low between 249–240 cm, 202–104 cm, and 75–64 cm. No carbonate was detected by LOI analysis and thus not reported here. Molybdenum (Mo) values range from 0–23 ppm with an average of 8.1 ± 6.7 ppm (Fig. 2). Mo is virtually absent or low between 249–246 cm and 190–135 cm. There is a notable decrease in Mo between 75–65 cm. Titanium (Ti) values range from 0–3384 with an average of 1184.8 ± 1117 ppm (Fig. 2). Ti is absent or very low between 260–251 cm, 221–212 cm, 107–79 cm, and 57–48 cm. Percent clay ranges from 4.5–21.6% with an average of 14.2 ± 3.5% (Fig. 2). Clay is highly variable between 260–166 cm, uniformly high between 166–132 cm, low and variable between 125–79 cm, and above average between 79–48 cm. Percent silt ranges from 59.4–87.2 % with an average of 78.1 ± 4.0% (Fig. 2). Silt is highly variable throughout the core with above-average values between 132–94 cm. Finally, percent sand ranges from 0.4–32.4 % with an average of 7.6 ± 4.9% (Fig. 2). Sand is variable and generally above average between 260–165 cm and 97–84 cm. Table 2 shows the various correlation coefficients (r-values) and their significance (P-values) for the data above.

Table 2. Sediment property correlation coefficients (r-value) and significance (P-value).

*bold significant at P < 0.0001

**italics significant at P < 0.05

Statistical analysis of the dataset revealed five coherent groups/clusters at the P < 0.001 within the dataset (Fig. 4). The PCA plot depicts that samples with high molybdenum (Mo) concentrations and high percent total organic matter (%TOM) occupy negative PC1 axis values (left); whereas, positive PC1 axis values (right) are dominated by high magnetic susceptibility (MS), high titanium (Ti) concentrations, and somewhat higher percent clay. PC1 axis accounts for 66.7 % of the variation within the samples. In contrast, grain size (i.e., sand vs. clay) controls the variation along the PC2 axis (i.e., 20.8% of the variation). Samples higher in percent sand plot within the negative PC2 axis values, while samples higher in percent clay plot in the positive PC2 axis range. Together, PC1 and PC2 account for 87.5% of the sample variation. However, we focus on PC1 only as the predominant integrated signal for changes in the lake's relative depth over time (see Discussion). PC2 is not addressed because it accounts for only 20.8% of the variance, and it did not add significant (or meaningful) information to our interpretation.

Figure 4. Principal component analysis (PCA) scatter plot of the sample depths with measurements for MS, %TOM, Mo, Ti, %clay, and %sand. All samples (symbols) and variables (lines) are plotted with respect to the first two eigenvectors (PC1 and PC2) determined from the PCA. The symbols represent sample groups that differ significantly (P-value <0.001) based on the SIMPROF analysis. We coded the samples according to these five clusters on the PCA plot to aid in the interpretation of our data: C1 blue (inverted triangles); C2 orange (normal triangles); C3 red (diamonds, squares, and circles). The latter symbols (diamonds, squares, and circles) were grouped together and color coded red because they represent only 14 of the 211 samples analyzed (i.e., <7 % of the population) and predominantly reflect small changes in percent sand. The clusters were color coded to show relative lake level (i.e., blue = relatively deep water; orange = relatively shallow water; red = transitional or variable relative lake level).

DISCUSSION

Sediment interpretations and statistics

Lake depth and productivity indicators

Considering its nearshore location, the core analyzed in this study is interpreted to reflect primarily a relative lake level signal and its effect on the lake's sedimentology and productivity. Nearshore core sites are valuable for capturing lake level transgressions and regressions because the nearshore environment is where the largest sedimentological changes usually occur in response to lake level fluctuations (Lehman, 1975; Stanley and Wear, 1978; Håkanson and Jansson, 1983; Dearing, 1991; R.S. Anderson et al., 2008; Pribyl and Shuman, 2014; Kirby et al., 2015; Bird et al., 2017; Leidelmeijer et al., 2021). We acknowledge that nearshore core sites are also subject to potential hiatuses and/or reductions in sediment preservation caused by regression-related erosion. In fact, we observed a significant reduction in sedimentation rates—possibly a hiatus—associated with a long-term Middle to Late Holocene regression. Despite the caveats associated with cores from a nearshore environment, the sedimentological, elemental, and biological data from Maddox Lake reveal large changes that are congruent with our relative lake level interpretation (see below). We also note that we use the term ‘relative lake level’ because we cannot assign absolute depths to our indicator-based, lake level history. Moreover, our shallow water core site, including the 1-m above modern surface elevation outlet, dictates that water depth likely did not exceed 3.6 m total depth at our core location (i.e., 2.6 m core plus 1 m outlet elevation). As a result, we use the terms ‘relative deep water’ or highstand and ‘relative shallow water’ or lowstand knowing that maximum water depth at our core site was never more than 3.6 m depth and declined in relative maximum depth at the core site as the basin's accommodation space diminished. Importantly, water depth need not change much (<<1 m in some cases) to impart a distinct and measurable sedimentological signature. For example, studies from wetlands, bogs, shallow lakes, marshes, littoral zone lake cores, etc., often contain rich sedimentological histories imparted by small changes in water depth (e.g., Mullins, 1998; Kirby et al., 2002; Mensing et al., 2013; Pigati et al., 2014; Wahl et al., 2016; Leeper et al., 2017; Honke et al., 2019; Longman et al., 2019; Anderson et al., 2022).

Changes in lake sediment grain size are often related to lake depth/lake level, distance from the shoreline, and grain size/energy dynamics associated with runoff from the drainage basin (Lehman, 1975; Anderson, 1977; Stanley and Wear, 1978; Håkanson and Jansson, 1983; Dearing, 1991; Anderson et al., 2008; Kirby et al., 2010, 2015, 2018; Pribyl and Shuman, 2014; Bird et al., 2017; Shuman and Serravezza, 2017; Leidelmeijer et al., 2021). In general, finer grain sizes are deposited in lower energy environments, such as a lake's depocenter or its most distal point from the shoreline. Conversely, coarser sediments are associated with the higher energy shoreline environments, where the finer sediment is easily eroded and focused into the deeper water (Lehman, 1975; Håkanson and Jansson, 1983; Dearing, 1991; Kirby et al., 2015; Bird et al., 2017; Pribyl and Shuman, 2014).

At Maddox Lake, clay and sand are negatively correlated (r = −0.59, P < 0.0001), as would be expected if changes in lake level also produce a change in the nearshore energy dynamics associated with the deposition of clay (deeper water = lower energy) versus sand (shallower water = higher energy) (Håkanson and Jansson, 1983; Pribyl and Shuman, 2014; Leidelmeijer et al., 2021) (Table 2). Periods of higher relative lake level also should be associated with more runoff into the basin. This is especially relevant where changes in winter precipitation and its associated runoff are the predominant controls on a lake's hydrology (Kirby et al., 2010, 2012, 2014; Hiner et al., 2016). Summer evaporation is also important in the context of changing Holocene summer insolation. Commonly used proxies for runoff are Ti and magnetic susceptibility (Thompson et al., 1975; Dearing, 1997; Haug et al., 2001; Brown et al., 2002; Kirby et al., 2004; Martin-Puertas et al., 2012). Positive correlations between Ti (ppm) and clay (r = 0.47, P < 0.0001) and MS and clay (r = 0.55, P < 0.0001) suggest that periods of lower energy at the core site (higher clay = deeper water) also are associated with higher Ti and MS (i.e., more runoff) (Table 2). Ti and MS are also positively correlated (r = 0.92, P < 0.0001), suggesting that both sediment components are recording the same depositional signal (i.e., runoff) (Table 2). Although not shown as plotted data, both Fe and Mn are also positively and significantly correlated with Ti and MS, suggesting a runoff and/or weathering source for Fe and Mn (Ti and Mn, r = 0.96, P < 0.0001; Ti and Fe, r = 0.99, P < 0.0001; MS and Fe, r = 0.93, P < 0.0001; MS and Mn, r = 0.93, P < 0.0001). In summary, we interpret higher clay content, Ti, Fe, Mn, and MS as indicators of a deeper relative lake with enhanced runoff (Table 3).

Table 3. Lake level interpretations.

Periods of lower lake level and less runoff are associated with lower clay content, Ti, and MS and with higher %TOM, oogonia counts, Mo, and sand content. For example, clay is negatively correlated with percent sand (r = −0.59, P < 0.0001), suggesting a response to nearshore energy dynamics. As lake level drops, coarser sediment is preferentially stored in the nearshore environment while the finer sediment is focused into the deeper basin. Clay is also negatively correlated with %TOM (r = −0.54, P < 0.0001) (Table 2). We interpret this relationship to reflect an increase in macrophyte growth and abundance in the nearshore environment during lower lake levels. As a modern analog, Maddox Lake is presently a shallow lake and macrophytes are abundant across the basin. In fact, modern surface sediments from the core location are characterized by high %TOM (78–80 %), high Mo (23–24 ppm), low Ti (0 ppm), and low MS (0 to −0.1). The positive relationship between %TOM and oogonia counts (r = 0.53, P < 0.0001) (Table 2) also supports this relative lake level interpretation. As lake level drops, Chara (sp.) proliferate in the nearshore environment increasing the abundance of oogonia and their preservation in the sediment column. The reduction or absence of oogonia during relatively deeper-water lake states may be the product of Chara (sp.) migration shoreward of the core site, enhanced turbidity in the photic zone due to more runoff, limited preservation in deeper water (which is potentially more caustic to calcitic preservation), or some combination of those factors (Spence, 1982; Andrews et al., 1984). The %TOM and oogonia counts are also negatively correlated to MS (r = −0.87, P < 0.0001; r = −0.55, P < 0.0001, respectively) (Table 2). Magnetic minerals are subject to dissolution under organic-rich, reducing conditions (Karlin and Levi, 1983; Hilton and Lishman, 1985; Canfield and Berner, 1987; Anderson and Rippey, 1988; Tarduno, 1995; Reynolds et al., 1999; Leidelmeijer et al., 2021). As the lake shoals and macrophytes thrive in the shallow water environment, magnetic minerals may experience selective dissolution. Alternatively, less runoff during lower lake levels (i.e., drier climates) may reduce the flux of magnetic minerals into the lake basin, as supported by low Ti, low Fe, and low Mn during inferred relative lowstands. In either case, lower MS values in the context of higher %TOM and oogonia counts are interpreted here to reflect lower relative lake levels. Finally, Mo is negatively correlated with MS (r = −0.84, P < 0.0001), clay content (r = −0.48, P < 0.0001), Fe (r = −0.91, P < 0.0001), and Mn (r = −0.89, P < 0.0001) and positively correlated with %TOM (r = 0.82, P < 0.0001) (Table 2). Mo is often used as an indicator for sediment column or hypolimnic anoxia and its concentration is often positively correlated to organic matter in lacustrine settings (Adelson et al., 2001; Eusterhues et al., 2005; Whitlock et al., 2012; Dahl et al., 2013). As lake level drops across the core site, the concentration of Mo should increase in the presence of higher organic matter (i.e., more oxygen consumption via decomposition/more reducing conditions) as the aquatic macrophytes dominate the nearshore environment—similar to the modern lake. The negative correlation of Mo with MS provides support for our interpretation for magnetic mineral dissolution during lowstands, when productivity is higher and sediment-water interface favors reducing conditions. However, the strong negative and statistically significant correlations between Mo-%TOM and Fe-Mn-Ti suggest that runoff is reduced during periods of inferred shallower water, when the lake is characterized by higher productivity. In either scenario, the data, when considered together, suggest that higher %TOM, oogonia counts, Mo, and sand, and lower Ti, Fe, Mn, MS, and clay content are evidence for a shallow, productive lake with less runoff (Table 3).

Statistics

To further assess these relationships with the objective to develop an integrated lake level indicator, we used multivariate statistical analyses. PC1 accounts for 66.7% of the variance (Fig. 4). This variance is caused by the opposing relationships between Mo and %TOM (to a lesser extent percent sand) and MS, Ti, and percent clay. Supported by simple correlation analyses, we conclude that PC1 best represents an integrated relative lake depth signal. Thus, we interpret positive PC1 values as a relatively deep-water lake environment (i.e., wetter climate) and negative PC1 values as a relatively shallow-water lake environment (i.e., drier climate) (Fig. 5). Most importantly, using this multivariate statistical methodology removes the subjective determination often associated with assigning units in a sediment core.

Figure 5. Core MLRC18-1 sediment data versus calibrated age. From bottom to top: (A) Magnetic susceptibility (× 10^-7 m³/kg). Core depth shown above the bottom x-axis, (B) number of oogonia capsules per 1 g dry sediment >125 μm, (C) percent total organic matter, (D) Mo concentration (ppm), (E) Ti concentration (ppm), (F) percent clay, (G) percent silt, (H) percent sand. C1 = Cluster 1 (blue), C2 = Cluster 2 (orange), C3 = Cluster 3 (red) (see Fig 4). LIA = Little Ice Age, MCA = Medieval Climatic Anomaly, ISI = Ice sheet influence.

Using our statistically determined lake level model, we plot the sediment data, including PC1, on Figure 5. The data are divided into three distinct clusters (C1, C2, and C3) based on our statistical analysis. C1 (blue) represents the deep lake units (i.e., high MS, high Ti, above average clay), C2 (orange) represents the shallow lake units (i.e., high TOM, high Mo) and C3 (red) represents the transitional lake units (variable sand). Notably, C3 represents only 14 of the 211 samples analyzed (i.e., <7% of the population) and predominantly reflects small changes in percent sand and clay. To avoid over-interpretation of these 14 data points, we do not develop C3 interpretatively, except to suggest these data represent transition stages in the lake's history. Plotting these clusters against the raw data and PC1 reveals large changes in the lake's sedimentology/limnology and relative lake level over time.

Maddox Lake relative 9000-year lake level history and climatic forcings

Variable Early-to-Middle Holocene lake level and the 8.2 ka cold event relative highstand

The Early Holocene at Maddox Lake is characterized by variably low lake levels from 9.0–7.4 cal ka BP, bracketing a pronounced highstand between 8.4–8.06 cal ka BP. We attribute this variably dry Early Holocene to the persistence of the Cordilleran and Laurentide Ice Sheets and its influence on the steering of winter storm tracks south of the study site. Additional forcings may include higher summer insolation (i.e., greater summer evaporation), a weaker Early Holocene latitudinal thermal gradient (i.e., weaker westerlies = less vigorous/less frequent winter storms across the study site), and warmer northeast Pacific SSTs (i.e., positive PDO-like conditions = less winter precipitation across the study site) (Barron et al., 2003; Carlson et al., 2008; Barron and Anderson, 2011; Steponaitis et al., 2015; Routson et al., 2019) (Fig. 6). Together, these Early Holocene forcings acted to decrease the frequency of winter storms across the study region, enhance summer evaporation, and sustain variably low lake level conditions at Maddox Lake.

Figure 6. Core MLRC18-1 PCA 1 relative level depth versus forcing data. (A) All left axes = PCA 1, (B) Difference in summer–winter insolation at 40 N latitude (w/m²) (Laskar et al., 1993), (C) ODP Site 1019 sea surface temperatures (°C) (Barron et al., 2003), (D) Cariaco Basin percent Ti (Haug et al., 2001), (E) NGRIP δ¹⁸O_(ice) (‰) (Rasmussen et al., 2006), (F) Western Pacific warm pool (WPWP)–Eastern Pacific warm pool (EPWP) tropical Pacific SST gradient anomaly (Koutavas and Joanides, 2012). C1 = Cluster 1 (blue), C2 = Cluster 2 (orange), C3 = Cluster 3 (red) (see figure 4). LIA = Little Ice Age, MCA = Medieval Climatic Anomaly, ISI = Ice sheet influence. Figure A/C and A/F PCA 1 show the 95 % range for the Bacon age model output in years max. and min. years from the median age.

Superimposed on this Early Holocene period of variably low lake levels is an abrupt and short-lived relative highstand between 8.4–8.06 cal ka BP (Fig. 5), which includes the 8.2 ka cold event (Alley et al., 1997; Barber et al., 1999; Lewis et al., 2012). Although recognized throughout the Northern Hemisphere, there are no records reporting this event from northwest CA. The geographically closest records capturing the 8.2 ka cold event in CA are located at White Moon Cave (~400 km south of Maddox Lake; Oster et al., 2017; de Wet et al., 2021) and carbonate deposits from Henry Cowell Park (~400 km south of Maddox Lake; Kanner et al., 2022) (Fig. 1). Oster et al. (2017) found evidence for more vigorous and/or more frequent winter storms across the region during the 8.2 ka cold event. From the same area, Kanner et al. (2022) found evidence from a perched tufa for wetter conditions ca. 8.0 ± 0.04 (2σ) cal ka BP. Our evidence for a lake level transgression (i.e., wetter winter climate) ca. the 8.2 ka cold event chronozone suggests a period of enhanced winter precipitation, which agrees with Oster et al. (2017), de Wet et al. (2021), and Kanner et al. (2022). Constrained by a lower calibrated ¹⁴C age (8.77 cal ka BP [8.81–8.64 cal ka BP 2σ]) and an upper calibrated ¹⁴C age (7.9 cal ka BP) [7.90–7.86 cal ka BP 2σ]), the Maddox Lake chronology (Bacon age model output = 8.4–8.06 cal ka BP) for the 8.2 ka event cannot definitively align our lake level transgression directly to the 8.2 ka cold event as defined in ice core chronologies (8.25–8.09 cal ka BP [Thomas et al., 2007]). However, other studies (Rohling and Palike, 2005) suggest that the 8.2 ka cold event was preceded by climatic instability beginning as early as 8.6 cal ka BP. As a result, the best we can say (based on our age model-constrained Maddox Lake transgression between 8.4–8.06 cal ka BP) is that this regression may correspond to the 8.2 ka cold event (see Fig. 2). Additional age control (beyond the scope of this paper) is required to confirm the regression's absolute timing relative to the 8.2 ka cold event, sensu stricto.

The 8.2 ka cold event generally is attributed to an abrupt increase in meltwater into the North Atlantic and the associated slowdown of the Atlantic Meridional Overturning Circulation (Barber et al., 1999; Renssen et al., 2001; Wiersma and Renssen, 2006; Morrill et al., 2013). Although 8.2 ka cold event climate models do not show significant changes to either temperature or precipitation along the CA coast (Wiersma and Renssen, 2006; Morrill et al., 2013), the 8.2 ka cold event is well documented throughout the Northern Hemisphere, suggesting that its effect was rapidly disseminated throughout the ocean-atmosphere system (e.g., Morrill and Jacobsen, 2005; Lutz et al., 2007; Nicolussi and Schlüchter, 2012; Young et al., 2012; Chabangborn et al., 2020). Climate models indicate a strengthening of the Aleutian Low and a cooling of the North Pacific in response to enhanced meltwater discharge into the North Atlantic (Renssen et al., 2001; Barron et al., 2003; Okumura et al., 2009; Wong et al., 2016). For example, the Maddox Lake 8.4–8.06 cal ka BP transgression is coeval with an abrupt and short-lived decrease in SSTs at the North Pacific ODP Site 1019 (Fig. 6). Together, these changes likely invigorated the westerlies and the subsequent frequency of winter season storm tracks across the study region (Okumura et al., 2009; Wong et al., 2016). Once the AMOC-driven 8.2 ka cold event forcings ceased, the lake returned to its pre-8.2 ka cold event relative lowstand.

Middle-to-Late Holocene regression and MCA relative lowstand

Following the variably dry Early Holocene, transitional lake states (C3) dominated the period 7.6–7.3 cal ka BP before returning to a relatively deep lake environment at 7.3 cal ka BP. The latter highstand lasted from 7.3 cal ka BP to ca. 5.2 cal ka BP. By 8.0–6.8 cal ka BP, the influence of the Cordilleran and Laurentide ice sheets as well as other glacial boundary conditions (e.g., greenhouses gas concentrations) was significantly diminished (Carlson et al., 2008; Routson et al., 2019). With the cessation of the Cordilleran and Laurentide Ice Sheet influence, other forcings such as insolation and Pacific ocean-atmosphere dynamics likely became more important.

A return to wetter conditions from 7.3–5.2 cal ka BP following the Early Holocene lowstand suggests more frequent winter storms across the study region. From 7.3–5.2 cal ka BP, the ITCZ was positioned farther north than modern, tropical Pacific SSTs favored La Niña conditions, and northeast Pacific SSTs were lower than today (negative PDO) (Fig. 6) (Haug et al., 2001; Barron et al., 2003; Koutavas and Joanides, 2012). Under modern climatological conditions, these ocean-atmosphere dynamics are generally associated with wetter winters north of the western US precipitation dipole (40°N [i.e., Maddox Lake location]) and drier conditions south (Wise, 2010, 2016). Summer insolation was also higher than modern from 7.3–5.2 cal ka BP, stressing the hydrologic budget of lake systems in the western US via enhanced summer evaporation (Steponaitis et al., 2015; Lachniet et al., 2020). This combination of more-frequent winter storms with enhanced summer evaporation may explain both the highstand's initiation at 7.3 cal ka BP and its transition by 5.2 cal ka BP to the Middle Holocene lake level regression. Interestingly, the hydroclimatic transition at ca. 5.2 cal ka BP is observed elsewhere in indicator-based reconstructions from North America (Forman et al., 2001; Shuman and Marsicek, 2016), Eurasia (Constantin et al., 2007; Bird et al., 2014; Robles et al., 2022), northwest Europe (Roland et al., 2015), South America (Thompson et al., 2006), and in Northern Hemisphere climate model simulations (Yahui et al., 2019).

The 7.3–5.2 cal ka BP relative highstand was followed by a subsequent lake level regression beginning by 5.2 cal ka BP. The lake became consistently shallow by 2.7 cal ka BP (Fig. 5). A persistent steady decline in PC1 values from 2.7–0.55 cal ka BP suggests a continued regression, culminating in the MCA (0.86–0.55 cal ka BP—according to the Maddox Lake age model). Notably, a dramatic decrease in sedimentation rates beginning at 6.8 cal ka BP suggests that the Middle Holocene regression begins earlier than 5.2 cal ka BP, when the lake was still deep (Fig. 2). Approximately 0.6 m of sediment was deposited between 6.8–0.86 cal ka BP. Although this interval is well dated (Table 1), we cannot rule out a potential hiatus sometime during this interval (6.8–0.86 cal ka BP, but likely after 5.2 cal ka BP based on the sedimentology) because the lake may have temporarily desiccated during the Late Holocene or MCA. There is no obvious stratigraphic or sedimentological evidence for desiccation; however, the low sedimentation rates and tightly spaced dates are curious and may reflect the absence of sediment via lowstand erosion rather than a decrease in sediment accumulation.

The long-term Middle-to-Late Holocene regression (beginning ca. 5.2 cal ka BP) likely records two signals: an internal geomorphic signal associated with the volumetric loss of basin accommodation space and an external climatic forcing signal. Absolute attribution to either the internal or external drivers cannot be resolved with the data available. Nonetheless, it is worth noting this interplay of forcings because of our nearshore core location and its hypersensitivity to changes in sedimentation. Over time, all lakes (excepting some active tectonic lakes) undergo a loss of accommodation space due to sedimentation. As previously noted, our nearshore core site likely did not exceed 3.6 m total depth (i.e., 2.6 m core plus 1-m outlet elevation) and declined in relative maximum depth as the basin's accommodation space diminished. This internal, geomorphic shoaling of the core site likely imparts some of the regressive sediment signal observed between 5.2–0.55 cal ka BP. However, measurable and non-linear changes in sediment properties within this period of shoaling also suggest non-geomorphic controls. Over the Middle-to-Late Holocene, the ITCZ migrated south, tropical Pacific SSTs transitioned to more frequent El Niño conditions, and northeast Pacific SSTs increased (positive PDO) (Fig. 6) (Haug et al., 2001; Barron et al., 2003; Koutavas and Joanides, 2012). This scenario often produces less winter precipitation north of the western US precipitation dipole (40°N) and wetter conditions south (Wise, 2010, 2016). Moreover, as Middle Holocene winter insolation decreased, the winter season latitudinal thermal gradient weakened, and winter storms may have tracked less frequently across the study region (Routson et al., 2019). Overall, our evidence for a persistent, lake level regression (i.e., climatic drying) at Maddox Lake over the Middle-to-Late Holocene suggests a western precipitation dipole that was positioned similar to modern or slightly southward and a similar-to-modern response to Pacific ocean-atmosphere dynamics.

LIA relative highstand and modern shallow lake

Superimposed on the millennial-scale Middle-to-Late Holocene Maddox Lake regression is a prominent Late Holocene highstand (0.54–0.37 cal ka BP), which we attribute to the LIA (Fig. 6) (Robock, 1979; Mann et al., 2009; Masson-Delmotte et al., 2013). Because this brief transgression is an obvious departure from the long-term lake level regression, we postulate that it reflects a change in the dominant climatic forcings at the time and cannot be explained via an internal geomorphic forcing. Although a matter of debate, North Atlantic volcanism and solar variability are the most common explanations for initiating the LIA (Mann et al., 2009; Miller et al., 2012; Slawinska and Robock, 2018). Climate models and indicator data suggest that the ocean-atmosphere coupled system rapidly responded to North Atlantic volcanism and solar variability before/during the LIA, resulting in a brief southward departure of the ITCZ, a possible reduction in AMOC, more El Niño-like conditions in the tropical Pacific, and possible expansion of Arctic sea ice (Haug et al., 2001; Chiang and Friedman, 2012; Koutavas and Joanides, 2012; Rustic et al., 2015; Cvijanovic et al., 2017; Slawinska and Robock, 2018; Lapointe and Bradley, 2021).

Throughout the Middle-to-Late Holocene, the ITCZ migrated south, north Pacific SSTs warmed, and El Niño-like conditions developed in the tropical Pacific (Fig. 6); yet Maddox Lake levels regressed in response to inferred less winter precipitation. How then, can we explain wetter conditions during the LIA under the same set of ocean-atmosphere dynamics? One possible explanation is that the LIA involved additional forcings such as a proposed increase in North Atlantic volcanism, solar variability, a reduction in AMOC, and expansion of Arctic sea ice (Mann et al., 2009; Chiang and Friedman, 2012; Miller et al., 2012; Cvijanovic et al., 2017; Slawinska and Robock, 2018; Lapointe and Bradley, 2021). These changes in volcanism, solar variability, and AMOC, coupled with an already southward positioned ITCZ and more El Niño-like conditions in the tropical Pacific, may have generated the dynamical “boost” required to shift the position of the western United States precipitation dipole northward. The result was a temporary increase in winter precipitation at Maddox Lake under conditions that normally favored drier winters. Clearly, the hydroclimatic system is complex and its response to far-field forcings spatially diverse. Additional records and regional modeling are required to analyze the finer-scale response of California precipitation to North Atlantic forcings throughout the Holocene, especially events such as the LIA that represent a departure from the predominant Holocene hydroclimatic trends.

This brief LIA relative highstand was followed by a return to low lake levels from 0.36–0.30 cal ka BP. Although not presented, the topmost section of the core (upper 48 cm) suggests continued low lake levels through to the modern as evinced by a thick peat unit with abundant, dense macrophytic roots. Notably, sedimentation rates were unusually high following the LIA highstand, which is counterintuitive to the proposed post-LIA lowstand interpretation. However, this peaty unit is uncompacted, water- and organic-rich, and lacking abundant clastic material. Thus, the sedimentation rate is exaggerated without accounting for density differences. In summary, the new Maddox Lake, relative lake level record captures key hydroclimatic changes recognized elsewhere in northwest CA. Briefly, we address the key hydroclimatic changes below in the context of the most-proximal records.

Some regional comparisons

Pollen-based vegetation reconstructions from northwest CA suggest that the Early Holocene was warmer and drier than present, the Middle Holocene was cooler and wetter than the Early Holocene, and the Late Holocene was drier than the Middle Holocene and drier than present (see Briles, 2017; Fig. 1). This general picture of vegetation-interpreted moisture approximates the same long-term decline in available moisture recorded at Maddox Lake. Unfortunately, there are no lake sediment studies similar to Maddox Lake that cover the past 9000 years from northwest CA for a more direct indicator-to-indicator comparison.

We also compared our relative lake level reconstruction to the 8000-year Oregon Caves National Monument (OCNM) δ¹⁸O_{(speleothem calcite)} atmospheric temperature record (Ersek et al., 2012) (Fig. 1). Located 172 km north of Maddox Lake, the OCNM represents a key hydroclimatic archive for regional comparisons in northern CA and southern OR (Fig. 7). In general, the OCNM δ¹⁸O_{(speleothem calcite)} data suggest a long-term increase in temperature from the Early to Late Holocene, although with substantial variability. This long-term increase in temperature may reflect a long-term increase in evaporation over the Holocene, fitting with the long-term lake level regression inferred at Maddox Lake. Overall, the lack of strong similarity between the OCNM and Maddox Lake records is not surprising given the very different proxies used, their disparate sample resolutions, and their varied response functions and sensitives.

Figure 7. Core MLRC18-1 PCA 1 relative level depth versus OCNM δ¹⁸O_{(speleothem calcite)} data (Ersek et al., 2012). C1 = Cluster 1 (blue), C2 = Cluster 2 (orange), C3 = Cluster 3 (red) (see Fig. 4). LIA = Little Ice Age, MCA = Medieval Climatic Anomaly, ISI = Ice sheet influence.

The Late Holocene Dry Period (LHDP: 2.8–1.85 cal ka BP) proposed by Mensing et al. (2013) represents an unusually pervasive period of dry climate across the Great Basin and the American West (Mensing et al., 2023). Although, continued research and retro-analysis of previously published papers and their age models indicate that the LHDP is of various durations and spatially diverse (e.g., Kirby et al., 2010, 2012, 2014, in southern CA). In northwest CA, however, the LHDP is not well documented. Even for this study (Maddox Lake), the LHDP is not a distinct period and more the culmination of a long-term decrease in winter precipitation and subsequent decline in lake level (Fig. 5).

Unlike the LHDP, the MCA (0.86–0.55 cal ka BP) at Maddox Lake is a distinct departure to low lake levels—perhaps the lowest lake level over the 9000-year record (Fig. 5). Elsewhere in northwest CA/southern OR, the MCA was characterized by large fires in the Upper Squaw Lake (OR) drainage basin (Colombaroli and Gavin, 2010) (Fig. 1). Mumbo Lake (CA), Bluff Lake (CA), and Crater Lake (CA) also show a period of higher-than-average fire activity between 0.8–1.2 cal ka BP (Mohr et al., 2000; Daniels et al., 2005) (Fig. 1). Data from Flycatcher Basin on the Modoc Plateau (CA) indicate a drier MCA based on a mesic versus xeric pollen ratio (Anderson et al., 2008) (Fig. 1). δ¹³C_{(speleothem calcite)} data suggest drier conditions during the MCA at OCNM (Ersek et al., 2012). A tree ring reconstruction of the Sacramento River flow was variable, but lower during the MCA (Meko et al., 2001).

The final significant change recorded at Maddox Lake occurred during the Little Ice Age—a brief but clear return to wet conditions (0.54–0.37 cal ka BP). Vegetation changes in the last 300 years at Sanger Lake (CA) and Bolan Lake (CA) suggest cooler (wetter?) conditions (Briles et al., 2008) (Fig. 1). A mesic versus xeric pollen ratio from Flycatcher Basin on the Modoc Plateau (CA) indicates wetter conditions during the LIA (Anderson et al., 2008) (Fig. 1). Speleothem δ¹³C_{(speleothem calcite)} data from OCNM also suggest a wet LIA (Ersek et al., 2012) (Fig. 1). Sacramento River flow based on tree ring reconstructions was higher during LIA (Meko et al., 2001).

CONCLUSIONS

We presented a new lake sediment record (Maddox Lake) from the northern California Coast Range. We used a multi-indicator approach coupled with 16 ¹⁴C dates to infer changes in relative lake level over the past 9000 years. In general, the Early Holocene was characterized by variably low lake levels with a brief, but pronounced, relative highstand (8.4–8.06 cal ka BP) possibly related to the 8.2 ka cold event. CA speleothem and tufa records show evidence for more-frequent and/or more-vigorous winter storms and fluvial activity at the same time as the brief, Early Holocene Maddox Lake relative highstand. Together, these records reflect a far-field response to North Atlantic forcing associated with enhanced meltwater into the formative region of the AMOC. A similar response (i.e., wetter climate) is observed during the LIA relative highstand at Maddox Lake, although volcanism and solar variability are the most likely drivers of the LIA. Wetter CA winters during both the 8.2 ka cold event and the LIA suggest a coupled North Atlantic–northwest California teleconnection, likely communicated through the coupled ocean-atmosphere system. Following a 7.3–5.2 cal ka BP highstand, lake levels declined consistently until 0.55 cal ka BP, culminating in the MCA relative lowstand. This gradual decrease in winter moisture availability during the Middle to Late Holocene is attributed to the combined effects of internal geomorphic forcings associated with loss of accommodation space and external forcings such as insolation, tropical and northeast Pacific SSTs, and the southward migration of the ITCZ. Together, these internal and external drivers acted to diminish the maximum water depth at the core site via sedimentation and volumetric loss and decrease the frequency/magnitude of winter storms tracking across the study region. Our new record provides evidence for a coupling between distal regions during the Holocene when the prevailing climate state was not too dissimilar from the modern. If the past is any predictor for the future, climate change and its effect on freshwater flux into the North Atlantic is likely to modulate CA's winter hydroclimatology, perhaps increasing the frequency and/or magnitude of winter storms.