INTRODUCTION
California faces a perennial freshwater availability crisis (MacDonald et al., Reference MacDonald, Kremenetski and Hidalgo2008; Diffenbaugh et al., Reference Diffenbaugh, Swain and Touma2015; Hatchett and McEvoy, Reference Hatchett and McEvoy2018). Key to the crisis is the feast-or-famine, winter-dominated climatology (Wang et al., Reference Wang, Yoon, Becker and Gillies2017; Swain et al., Reference Swain, Langenbrunner, Neelin and Hall2018). 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., Reference Das, Maurer, Pierce, Dettinger and Cayan2013; Swain et al., Reference Swain, Langenbrunner, Neelin and Hall2018; Ullrich et al., Reference Ullrich, Xu, Rhoades, Dettinger, Mount, Jones and Vahmani2018). 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., Reference Schmidt, Shindell, Miller, Mann and Rind2004; Shuman and Marsicek, Reference Shuman and Marsicek2016; Bakker et al., Reference Bakker, Clark, Golledge, Schmittner and Weber2017; Cook et al., Reference Cook, Mankin and Anchukaitis2018; Lachniet et al., Reference Lachniet, Asmerom, Polyak and Denniston2020; Tierney et al., Reference Tierney, Poulsen, Montañez, Bhattacharya, Feng, Ford, Hönisch, Inglis, Petersen and Sagoo2020). 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., Reference Barron, Heusser, Herbert and Lyle2003), 2. TN062-O550 (Barron et al., Reference Barron, Bukry, Heusser, Addison and Alexander2018), 3. Oregon Caves National Monument (Ersek et al., Reference Ersek, Clark, Mix, Cheng and Edwards2012), 4. Upper Squaw Lake (OR) (Colombaroli and Gavin, Reference Colombaroli and Gavin2010), 5. Sanger Lake (Briles et al., Reference Briles, Whitlock, Bartlein and Higuera2008, Reference Briles, Whitlock, Skinner and Mohr2011; Briles, Reference Briles2017), 6. Bolan Lake (OR) (Briles et al., Reference Briles, Whitlock and Bartlein2005, Reference Briles, Whitlock, Bartlein and Higuera2008; Whitlock et al., Reference Whitlock, Marlon, Briles, Brunelle, Long and Bartlein2008), 7. Twin Lakes (Wanket, Reference Wanket2002) and Fish Lake (Crawford et al., Reference Crawford, Mensing, Lake and Zimmerman2015), 8. Lake Ogaromtoc (Crawford et al., Reference Crawford, Mensing, Lake and Zimmerman2015), 9. Campbell Lake (Briles et al., Reference Briles, Whitlock, Skinner and Mohr2011; Briles, Reference Briles2017) and Taylor Lake (Briles et al., Reference Briles, Whitlock, Skinner and Mohr2011; Briles, Reference Briles2017), 10. Crater Lake (CA) (Mohr et al., Reference Mohr, Whitlock and Skinner2000), Bluff Lake (Mohr et al., Reference Mohr, Whitlock and Skinner2000), Cedar Lake (Briles et al., Reference Briles, Whitlock, Skinner and Mohr2011; Briles, Reference Briles2017), and Mumbo Lake (Daniels et al., Reference Daniels, Anderson and Whitlock2005), 11. Flycatcher Basin (R.S. Anderson et al., Reference Anderson, Smith, Jass and Spaulding2008), 12. White Moon Cave (Oster et al., Reference Oster, Sharp, Covey, Gibson, Rogers and Mix2017).
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, Reference Cayan, Peterson and Peterson1989; Cayan et al., Reference Cayan, Dettinger, Diaz and Graham1998; Dettinger et al., Reference Dettinger, Cayan, Diaz and Meko1998; Wise, Reference Wise2010). 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., Reference Ralph, Neiman, Wick, Gutman, Dettinger, Cayan and White2006; Dettinger et al., Reference Dettinger, Ralph, Das, Neiman and Cayan2011).
Superimposed on this unimodal winter climatology is the western United States precipitation dipole (Cayan et al., Reference Cayan, Dettinger, Diaz and Graham1998; Dettinger et al., Reference Dettinger, Cayan, Diaz and Meko1998; Fye et al., Reference Fye, Stahle and Cook2004; Wise, Reference Wise2010, Reference Wise2016). 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, Reference Wise2010, Reference Wise2016). Of course, this boundary is dynamic and occasionally absent altogether (Cayan et al., Reference Cayan, Dettinger, Diaz and Graham1998; Dettinger et al., Reference Dettinger, Cayan, Diaz and Meko1998; Fye et al., Reference Fye, Stahle and Cook2004; Wise, Reference Wise2010, Reference Wise2016). 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., Reference Benson, Lund, Paillet, Smoot, Kester, Mensing, Meko, Lindström, Kashgarian and Rye2002; Kirby et al., Reference Kirby, Poulsen, Lund, Patterson, Reidy and Hammond2004, Reference Kirby, Lund and Bird2006; Hanson et al., Reference Hanson, Dettinger and Newhouse2006; Adams et al., Reference Adams, Negrini, Cook and Rajagopal2015).
Wise (Reference Wise2010) 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, Reference Wise2010). 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., Reference Cayan, Dettinger, Diaz and Graham1998, Reference Cayan, Redmond and Riddle1999; Dettinger et al., Reference Dettinger, Cayan, Diaz and Meko1998; Wise, Reference Wise2010; Peng et al., Reference Peng, Yu and Gautam2013). 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, Reference Seager and Hoerling2014). 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., Reference Seager, Harnik, Robinson, Kushnir, Ting, Huang and Velez2005, Reference Seager, Hoerling, Schubert, Wang, Lyon, Kumar, Nakamura and Henderson2014; Cook et al., Reference Cook, Seager and Miller2011, Reference Cook, Seager, Miller and Mason2013; Chiang and Friedman, Reference Chiang and Friedman2012; Peng et al., Reference Peng, Yu and Gautam2013; Coats et al., Reference Coats, Smerdon, Cook, Seager, Cook and Anchukaitis2016; Kam and Sheffield, Reference Kam and Sheffield2016; Wise, Reference Wise2016). In the Atlantic, the AMO is an important player in modulating North Atlantic climate. As Wise (Reference Wise2010) and Zhang and Delworth (Reference Zhang and Delworth2007) demonstrated, the AMO affects Pacific climate as well. Wise (Reference Wise2010) 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 km2 with a drainage basin of ~1.54 km2. 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, Reference Adam and West1983; West, Reference West, White, Mikkelson, Hildebrandt, Hildebrandt and Basgall1993; Heusser, Reference Heusser1998; Mohr et al., Reference Mohr, Whitlock and Skinner2000; Meko et al., Reference Meko, Therrell, Baisan and Hughes2001; Barron et al., Reference Barron, Heusser, Herbert and Lyle2003; Daniels et al., Reference Daniels, Anderson and Whitlock2005; Vacco et al., Reference Vacco, Clark, Mix, Cheng and Edwards2005; Briles et al., Reference Briles, Whitlock, Bartlein and Higuera2008, Reference Briles, Whitlock, Skinner and Mohr2011; Whitlock et al., Reference Whitlock, Marlon, Briles, Brunelle, Long and Bartlein2008; Colombaroli and Gavin, Reference Colombaroli and Gavin2010; Ersek et al., Reference Ersek, Clark, Mix, Cheng and Edwards2012; Malevich et al., Reference Malevich, Woodhouse and Meko2013; Crawford et al., Reference Crawford, Mensing, Lake and Zimmerman2015; Briles, Reference Briles2017; Anderson et al., Reference Anderson, Wahl and Bhattacharya2020). 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, Reference Briles2017). Despite some of these complex responses, there emerges a general picture of Holocene fire and vegetation dynamics.
Briles (Reference Briles2017) 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, Reference Briles2017). A closer look at the Late Holocene from two of these seven study lakes (Fish Lake and Lake Ogaromtoc) by Crawford et al. (Reference Crawford, Mensing, Lake and Zimmerman2015) 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, Reference Colombaroli and Gavin2010).
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., Reference Ersek, Clark, Mix, Cheng and Edwards2012). Using stable isotopes δ13C(speleothem calcite) and δ18O(speleothem calcite), Ersek et al. (Reference Ersek, Clark, Mix, Cheng and Edwards2012) inferred winter precipitation amount and/or surface biomass development and atmospheric temperature, respectively. Low-frequency Holocene δ18O(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 δ13C(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. (Reference Ersek, Clark, Mix, Cheng and Edwards2012) attributed this higher frequency change to a combination of solar forcing and Pacific Ocean-atmosphere dynamics.
From the adjacent northeast Pacific, Barron et al. (Reference Barron, Heusser, Herbert and Lyle2003, Reference Barron, Bukry, Heusser, Addison and Alexander2018) 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., Reference Barron, Heusser, Herbert and Lyle2003). 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., Reference Barron, Heusser, Herbert and Lyle2003).
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., Reference Barron, Bukry, Heusser, Addison and Alexander2018). 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., Reference Barron, Heusser, Herbert and Lyle2003, Reference Barron, Bukry, Heusser, Addison and Alexander2018; Barron and Anderson, Reference Barron and Anderson2011).
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, Reference Harrison and Digerfeldt1993; Dearing, Reference Dearing1997; Shuman and Serravezza, Reference Shuman and Serravezza2017). 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] 14C 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 14C 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, D14C, and conventional radiocarbon age, following the conventions of Stuiver and Polach (Reference Stuiver and Polach1977, p. 355). Sample preparation backgrounds have been subtracted, based on measurements of 14C-free wood. All results have been corrected for isotopic fractionation according to the conventions of Stuiver and Polach (Reference Stuiver and Polach1977), with δ13C values measured on prepared graphite using the AMS spectrometer (these can differ from δ13C of the original material and are not shown). Samples labeled “Modern” contain excess 14C, 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 m3/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, Reference Dean1974). Mass accumulation rates were calculated by dividing the sedimentation rate per depth by the dry bulk density per depth (g/yr/cm2).
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., Reference Niederman, Porinchu and Kotlia2021; Wright et al., Reference Wright, Bird, Gibson, Pollard, Escobar and Barr2023). However, the data should be considered semi-quantitative due to inherent sediment heterogeneity and its variable clastic, organic, and chemical components (Weltje and Tjallingii, Reference Weltje and Tjallingii2008). 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, Reference Groves and Bullock-Webster1924; Burne et al., Reference Burne, Bauld and De Deckker1980; Vance et al., Reference Vance, Mathewes and Clague1992; Mullins, Reference Mullins1998; Apolinarska and Hammarlund, Reference Apolinarska and Hammarlund2009; Détriché et al., Reference Détriché, Bréhéret, Soulié-Märsche, Karrat and Macaire2009; Kemp et al., Reference Kemp, Radke, Olley, Juggins and De Deckker2012). 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% H2O2, 10 mL of 1N HCl, and 10 mL of 1N NaOH (Leidelmeijer et al., Reference Leidelmeijer, Kirby, MacDonald, Carlin, Avila, Han, Nauman, Loyd, Nichols and Ramezan2021). 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, Reference Wentworth1922) 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., Reference Clarke, Somerfield and Gorley2008, and Somerfield and Clarke, Reference Somerfield and Clarke2013). 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 14C ages (far left side; values given in cal yr BP). From left to right: (A) Magnetic susceptibility (× 10-7 m3/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” 14C 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, Reference Blaauw and Christen2011). 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, Reference Blaauw and Christen2011). Blue features are the calibrated 14C 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 m3/kg with an average of 0.4 ± 0.4 × 10-7 m3/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, Reference Lehman1975; Stanley and Wear, Reference Stanley and Wear1978; Håkanson and Jansson, Reference Håkanson and Jansson1983; Dearing, Reference Dearing1991; R.S. Anderson et al., Reference Anderson, Smith, Jass and Spaulding2008; Pribyl and Shuman, Reference Pribyl and Shuman2014; Kirby et al., Reference Kirby, Knell, Anderson, Lachniet, Palermo, Eeg and Lucero2015; Bird et al., Reference Bird, Lei, Perello, Polissar, Yao, Finney, Bain, Pompeani and Thompson2017; Leidelmeijer et al., Reference Leidelmeijer, Kirby, MacDonald, Carlin, Avila, Han, Nauman, Loyd, Nichols and Ramezan2021). 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, Reference Mullins1998; Kirby et al., Reference Kirby, Mullins, Patterson and Burnett2002; Mensing et al., Reference Mensing, Sharpe, Tunno, Sada, Thomas, Starratt and Smith2013; Pigati et al., Reference Pigati, Rech, Quade and Bright2014; Wahl et al., Reference Wahl, Hansen, Byrne, Anderson and Schreiner2016; Leeper et al., Reference Leeper, Rhodes, Kirby, Scharer, Carlin, Hemphill-Haley, Avnaim-Katav, MacDonald, Starratt and Aranda2017; Honke et al., Reference Honke, Pigati, Wilson, Bright, Goldstein, Skipp, Reheis and Havens2019; Longman et al., Reference Longman, Veres, Ersek, Haliuc and Wennrich2019; Anderson et al., Reference Anderson, Skipp, Strickland, Honke, Havens and VanSistine2022).
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, Reference Lehman1975; Anderson, Reference Anderson1977; Stanley and Wear, Reference Stanley and Wear1978; Håkanson and Jansson, Reference Håkanson and Jansson1983; Dearing, Reference Dearing1991; Anderson et al., Reference Anderson, Smith, Jass and Spaulding2008; Kirby et al., Reference Kirby, Lund, Patterson, Anderson, Bird, Ivanovici, Monarrez and Nielsen2010, Reference Kirby, Knell, Anderson, Lachniet, Palermo, Eeg and Lucero2015, Reference Kirby, Heusser, Scholz, Ramezan, Anderson, Markle and Rhodes2018; Pribyl and Shuman, Reference Pribyl and Shuman2014; Bird et al., Reference Bird, Lei, Perello, Polissar, Yao, Finney, Bain, Pompeani and Thompson2017; Shuman and Serravezza, Reference Shuman and Serravezza2017; Leidelmeijer et al., Reference Leidelmeijer, Kirby, MacDonald, Carlin, Avila, Han, Nauman, Loyd, Nichols and Ramezan2021). 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, Reference Lehman1975; Håkanson and Jansson, Reference Håkanson and Jansson1983; Dearing, Reference Dearing1991; Kirby et al., Reference Kirby, Knell, Anderson, Lachniet, Palermo, Eeg and Lucero2015; Bird et al., Reference Bird, Lei, Perello, Polissar, Yao, Finney, Bain, Pompeani and Thompson2017; Pribyl and Shuman, Reference Pribyl and Shuman2014).
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, Reference Håkanson and Jansson1983; Pribyl and Shuman, Reference Pribyl and Shuman2014; Leidelmeijer et al., Reference Leidelmeijer, Kirby, MacDonald, Carlin, Avila, Han, Nauman, Loyd, Nichols and Ramezan2021) (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., Reference Kirby, Lund, Patterson, Anderson, Bird, Ivanovici, Monarrez and Nielsen2010, Reference Kirby, Zimmerman, Patterson and Rivera2012, Reference Kirby, Feakins, Hiner, Fantozzi, Zimmerman, Dingemans and Mensing2014; Hiner et al., Reference Hiner, Kirby, Bonuso, Patterson, Palermo and Silveira2016). 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., Reference Thompson, Battarbee, O'Sullivan and Oldfield1975; Dearing, Reference Dearing1997; Haug et al., Reference Haug, Hughen, Sigman, Peterson and Rohl2001; Brown et al., Reference Brown, Bierman, Lini, Davis and Southon2002; Kirby et al., Reference Kirby, Poulsen, Lund, Patterson, Reidy and Hammond2004; Martin-Puertas et al., Reference Martin-Puertas, Brauer, Dulski and Brademann2012). 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, Reference Spence1982; Andrews et al., Reference Andrews, Davison, Andrews and Raven1984). 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, Reference Karlin and Levi1983; Hilton and Lishman, Reference Hilton and Lishman1985; Canfield and Berner, Reference Canfield and Berner1987; Anderson and Rippey, Reference Anderson and Rippey1988; Tarduno, Reference Tarduno1995; Reynolds et al., Reference Reynolds, Rosenbaum, van Metre, Tuttle, Callender and Goldin1999; Leidelmeijer et al., Reference Leidelmeijer, Kirby, MacDonald, Carlin, Avila, Han, Nauman, Loyd, Nichols and Ramezan2021). 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., Reference Adelson, Helz and Miller2001; Eusterhues et al., Reference Eusterhues, Heinrichs and Schneider2005; Whitlock et al., Reference Whitlock, Dean, Fritz, Stevens, Stone, Power, Rosenbaum, Pierce and Bracht-Flyr2012; Dahl et al., Reference Dahl, Ruhl, Hammarlund, Canfield, Rosing and Bjerrum2013). 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 m3/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., Reference Barron, Heusser, Herbert and Lyle2003; Carlson et al., Reference Carlson, Legrande, Oppo, Came, Schmidt, Anslow, Licciardi and Obbink2008; Barron and Anderson, Reference Barron and Anderson2011; Steponaitis et al., Reference Steponaitis, Andrews, McGee, Quade, Hsieh, Broecker, Shuman, Burns and Cheng2015; Routson et al., Reference Routson, McKay, Kaufman, Erb, Goosse, Shuman, Rodysill and Ault2019) (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/m2) (Laskar et al., Reference Laskar, Joutel and Boudin1993), (C) ODP Site 1019 sea surface temperatures (°C) (Barron et al., Reference Barron, Heusser, Herbert and Lyle2003), (D) Cariaco Basin percent Ti (Haug et al., Reference Haug, Hughen, Sigman, Peterson and Rohl2001), (E) NGRIP δ18O(ice) (‰) (Rasmussen et al., Reference Rasmussen, Andersen, Svensson, Steffensen, Vinther, Clausen and Siggaard-Andersen2006), (F) Western Pacific warm pool (WPWP)–Eastern Pacific warm pool (EPWP) tropical Pacific SST gradient anomaly (Koutavas and Joanides, Reference Koutavas and Joanides2012). 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., Reference Alley, Mayewski, Sowers, Stuiver, Taylor and Clark1997; Barber et al., Reference Barber, Dyke, Hillaire-Marcel, Jennings, Andrews, Kerwin and Bilodeau1999; Lewis et al., Reference Lewis, Miller, Levac, Piper and Sonnichsen2012). 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., Reference Oster, Sharp, Covey, Gibson, Rogers and Mix2017; de Wet et al., Reference de Wet, Erhardt, Sharp, Marks, Bradbury, Turchyn, Xu and Oster2021) and carbonate deposits from Henry Cowell Park (~400 km south of Maddox Lake; Kanner et al., Reference Kanner, Cortes and Ibarra2022) (Fig. 1). Oster et al. (Reference Oster, Sharp, Covey, Gibson, Rogers and Mix2017) 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. (Reference Kanner, Cortes and Ibarra2022) 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. (Reference Oster, Sharp, Covey, Gibson, Rogers and Mix2017), de Wet et al. (Reference de Wet, Erhardt, Sharp, Marks, Bradbury, Turchyn, Xu and Oster2021), and Kanner et al. (Reference Kanner, Cortes and Ibarra2022). Constrained by a lower calibrated 14C age (8.77 cal ka BP [8.81–8.64 cal ka BP 2σ]) and an upper calibrated 14C 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., Reference Thomas, Wolff, Mulvaney, Steffensen, Johnsen, Arrowsmith, White, Vaughn and Popp2007]). However, other studies (Rohling and Palike, Reference Rohling and Palike2005) 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., Reference Barber, Dyke, Hillaire-Marcel, Jennings, Andrews, Kerwin and Bilodeau1999; Renssen et al., Reference Renssen, Goosse, Fichefet and Campin2001; Wiersma and Renssen, Reference Wiersma and Renssen2006; Morrill et al., Reference Morrill, LeGrande, Renssen, Bakker and Otto-Bliesner2013). 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, Reference Wiersma and Renssen2006; Morrill et al., Reference Morrill, LeGrande, Renssen, Bakker and Otto-Bliesner2013), 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, Reference Morrill and Jacobsen2005; Lutz et al., Reference Lutz, Wiles, Lowell and Michaels2007; Nicolussi and Schlüchter, Reference Nicolussi and Schlüchter2012; Young et al., Reference Young, Briner, Rood and Finkel2012; Chabangborn et al., Reference Chabangborn, Punwong, Phountong, Nudnara, Yoojam, Sainakum, Won-In and Sompongchaiyakul2020). 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., Reference Renssen, Goosse, Fichefet and Campin2001; Barron et al., Reference Barron, Heusser, Herbert and Lyle2003; Okumura et al., Reference Okumura, Deser, Hu, Timmermann and Xie2009; Wong et al., Reference Wong, Potter, Montañez, Otto-Bliesner, Behling and Oster2016). 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., Reference Okumura, Deser, Hu, Timmermann and Xie2009; Wong et al., Reference Wong, Potter, Montañez, Otto-Bliesner, Behling and Oster2016). 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., Reference Carlson, Legrande, Oppo, Came, Schmidt, Anslow, Licciardi and Obbink2008; Routson et al., Reference Routson, McKay, Kaufman, Erb, Goosse, Shuman, Rodysill and Ault2019). 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., Reference Haug, Hughen, Sigman, Peterson and Rohl2001; Barron et al., Reference Barron, Heusser, Herbert and Lyle2003; Koutavas and Joanides, Reference Koutavas and Joanides2012). 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, Reference Wise2010, Reference Wise2016). 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., Reference Steponaitis, Andrews, McGee, Quade, Hsieh, Broecker, Shuman, Burns and Cheng2015; Lachniet et al., Reference Lachniet, Asmerom, Polyak and Denniston2020). 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., Reference Forman, Oglesby and Webb2001; Shuman and Marsicek, Reference Shuman and Marsicek2016), Eurasia (Constantin et al., Reference Constantin, Bojar, Lauritzen and Lundberg2007; Bird et al., Reference Bird, Polisar, Lei, Thompson, Yao, Finney, Bain, Pompeani and Steinman2014; Robles et al., Reference Robles, Peyron, Brugiapaglia, Ménot, Dugerdil, Ollivier and Ansanay-Alex2022), northwest Europe (Roland et al., Reference Roland, Daley, Caseldine, Charman, Turney, Amesbury, Thompson and Woodley2015), South America (Thompson et al., Reference Thompson, Mosley-Thompson, Brecher, Davis, León, Les, Lin, Mashiotta and Mountain2006), and in Northern Hemisphere climate model simulations (Yahui et al., Reference Yahui, Jian, Bin, Liang and Mi2019).
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., Reference Haug, Hughen, Sigman, Peterson and Rohl2001; Barron et al., Reference Barron, Heusser, Herbert and Lyle2003; Koutavas and Joanides, Reference Koutavas and Joanides2012). This scenario often produces less winter precipitation north of the western US precipitation dipole (40°N) and wetter conditions south (Wise, Reference Wise2010, Reference Wise2016). 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., Reference Routson, McKay, Kaufman, Erb, Goosse, Shuman, Rodysill and Ault2019). 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, Reference Robock1979; Mann et al., Reference Mann, Zhang, Rutherford, Bradley, Hughes, Shindell, Ammann, Faluvegi and Ni2009; Masson-Delmotte et al., Reference Masson-Delmotte, Schulz, Abe-Ouchi, Beer, Ganopolski, González Rouco, Jansen, Stocker, Qin, Plattner, Tignor, Allen, Boschung, Nauels, Xia, Bex and Midgley2013). 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., Reference Mann, Zhang, Rutherford, Bradley, Hughes, Shindell, Ammann, Faluvegi and Ni2009; Miller et al., Reference Miller, Geirsdottir, Zhong, Larsen, Otto-Bliesner, Holland and Bailey2012; Slawinska and Robock, Reference Slawinska and Robock2018). 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., Reference Haug, Hughen, Sigman, Peterson and Rohl2001; Chiang and Friedman, Reference Chiang and Friedman2012; Koutavas and Joanides, Reference Koutavas and Joanides2012; Rustic et al., Reference Rustic, Koutavas, Marchitto and Linsley2015; Cvijanovic et al., Reference Cvijanovic, Santer, Bonfils, Lucas, Chiang and Zimmerman2017; Slawinska and Robock, Reference Slawinska and Robock2018; Lapointe and Bradley, Reference Lapointe and Bradley2021).
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., Reference Mann, Zhang, Rutherford, Bradley, Hughes, Shindell, Ammann, Faluvegi and Ni2009; Chiang and Friedman, Reference Chiang and Friedman2012; Miller et al., Reference Miller, Geirsdottir, Zhong, Larsen, Otto-Bliesner, Holland and Bailey2012; Cvijanovic et al., Reference Cvijanovic, Santer, Bonfils, Lucas, Chiang and Zimmerman2017; Slawinska and Robock, Reference Slawinska and Robock2018; Lapointe and Bradley, Reference Lapointe and Bradley2021). 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, Reference Briles2017; 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) δ18O(speleothem calcite) atmospheric temperature record (Ersek et al., Reference Ersek, Clark, Mix, Cheng and Edwards2012) (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 δ18O(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 δ18O(speleothem calcite) data (Ersek et al., Reference Ersek, Clark, Mix, Cheng and Edwards2012). 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. (Reference Mensing, Sharpe, Tunno, Sada, Thomas, Starratt and Smith2013) represents an unusually pervasive period of dry climate across the Great Basin and the American West (Mensing et al., Reference Mensing, Wang, Rhode, Kennett, Csank, Thomas, Briem, Harper, Culleton and George2023). 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., Reference Kirby, Lund, Patterson, Anderson, Bird, Ivanovici, Monarrez and Nielsen2010, Reference Kirby, Zimmerman, Patterson and Rivera2012, Reference Kirby, Feakins, Hiner, Fantozzi, Zimmerman, Dingemans and Mensing2014, 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, Reference Colombaroli and Gavin2010) (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., Reference Mohr, Whitlock and Skinner2000; Daniels et al., Reference Daniels, Anderson and Whitlock2005) (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., Reference Anderson, Smith, Jass and Spaulding2008) (Fig. 1). δ13C(speleothem calcite) data suggest drier conditions during the MCA at OCNM (Ersek et al., Reference Ersek, Clark, Mix, Cheng and Edwards2012). A tree ring reconstruction of the Sacramento River flow was variable, but lower during the MCA (Meko et al., Reference Meko, Therrell, Baisan and Hughes2001).
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., Reference Briles, Whitlock, Bartlein and Higuera2008) (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., Reference Anderson, Smith, Jass and Spaulding2008) (Fig. 1). Speleothem δ13C(speleothem calcite) data from OCNM also suggest a wet LIA (Ersek et al., Reference Ersek, Clark, Mix, Cheng and Edwards2012) (Fig. 1). Sacramento River flow based on tree ring reconstructions was higher during LIA (Meko et al., Reference Meko, Therrell, Baisan and Hughes2001).
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 14C 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.
Acknowledgments
We acknowledge the native lands of the Northern Wintu, Nor Rel Muk Wintu, Cayuse, Umatilla, and Walla Walla indigenous people where Maddox Lake is located, and we thank them and their ancestors for accessing this resource (www.native-land.ca). Thanks also to the Shasta-Trinity National Forest and Dennis Veich (Forest Geologist), Lois Shoemaker, and Lesley Yen for providing site access and research permission. This research was funded by the National Science Foundation Grant to Kirby, Carlin, Nichols, Ramezan, and MacDonald (NSF-EAR #1702825). We thank the editor (Dr. Lesleigh Anderson) and two anonymous reviewers for their insightful and helpful edits and suggestions.
Open access
