During the Last Glacial Maximum (LGM) and subsequent deglaciation, the Great Basin in the southwestern United States was covered by numerous extensive closed-basin lakes, in stark contrast with the predominately arid climate observed today. This transition from lakes in the Late Pleistocene to modern aridity implies large changes in the regional water balance. Whether these changes were driven by increased precipitation rates due to changes in atmospheric dynamics, decreased evaporation rates resulting from temperature depression and summer insolation changes, or some combination of the two remains uncertain. The factors contributing to these large-scale changes in hydroclimate are critical to resolve, given that this region is poised to undergo future anthropogenic-forced climate changes with large uncertainties in model simulations for the 21st century. Furthermore, there are ambiguous constraints on the magnitude and even the sign of changes in key hydroclimate variables between the Last Glacial Maximum and the present day in both proxy reconstructions and climate model analyses of the region. Here we report thermodynamically derived estimates of changes in temperature, precipitation, and evaporation rates, as well as the isotopic composition of lake water, using clumped isotope data from an ancient lake in the northwestern Great Basin, Lake Surprise (California). Compared to modern climate, mean annual air temperature at Lake Surprise was 4.7 °C lower during the Last Glacial Maximum, with decreased evaporation rates and similar precipitation rates to modern. During the mid-deglacial period, the growth of Lake Surprise implied that the lake hydrologic budget briefly departed from steady state. Our reconstructions indicate that this growth took place rapidly, while the subsequent lake regression took place over several thousand years. Using models for precipitation and evaporation constrained from clumped isotope results, we determine that the disappearance of Lake Surprise coincided with a moderate increase in lake temperature, along with increasing evaporation rates outpacing increasing precipitation rates. Concomitant analysis of proxy data and climate model simulations for the Last Glacial Maximum are used to provide a robust means to understand past climate change, and by extension, predict how current hydroclimates may respond to expected future climate forcings. We suggest that an expansion of this analysis to more basins across a larger spatial scale could provide valuable insight into proposed climate forcings, and aid in climate model process depiction. Ultimately, our analysis highlights the importance of temperature-driven evaporation as a mechanism for lake growth and retreat in this region.


The Great Basin (southwestern United States) is characterized today by its aridity and low precipitation, with many regions receiving <250 mm of rain per year. In stark contrast, for the Last Glacial Maximum (LGM; ∼23,000–19,000 years ago; Clark et al., 2009) and subsequent deglaciation (∼19,000–11,000 years ago), the sedimentary record and landscape geomorphology indicate that the region was much wetter and marked by extensive lake systems in most endorheic basins (Ibarra et al., 2018; Mifflin and Wheat, 1979; Reheis, 1999). Where these ancient lakes once existed, dry salt flats (playas) now occupy the landscape. The abrupt change in regional hydroclimate implied by these playas is not unique to the LGM-to-Holocene transition or to the Great Basin. During recent decades, arid regions in many parts of the world have seen intensified drought, while many wet regions are reporting increased-intensity rainfall events (e.g., Rind et al., 1990; Zwiers and Kharin, 1998). On a global scale, the same period of time has also seen large-scale changes in zonal precipitation, whereby water-scarce latitudes are becoming drier, and wet latitudes are generally becoming wetter. It is commonly believed that this “rich get richer” mechanism is a result of global climate change and intensification of the hydrological cycle (e.g., Trenberth, 2011).

Anthropogenic climate change provides a substantial impetus to better understand the regional hydroclimate response to increasing air temperatures, because the stability of the current climate system has far-reaching implications that permeate through modern society. For example, with respect to water scarcity, the consequences of climate change can be severe—from economic hardships imposed on the agricultural industry, especially small-scale farming (e.g., Rosenweig and Hillel, 2008), to increasing wildfire occurrence (Abatzoglou and Williams, 2016) and disease prevalence (e.g., the West Nile virus; Epstein, 2001) to the preservation of native species (Belote et al., 2017)—and these effects are likely to be widespread.

One approach that can improve our understanding of different atmospheric processes that drive aridification in the Great Basin involves using paleoclimate data, in conjunction with data-model comparison, to study controls on past changes in the regional water balance. This approach is valuable not only because it allows for a better understanding of regional climate response to climate change, but also because it can shed light on the accuracy of climate model simulations of both past and future temperature, precipitation, and evaporation changes. Thus, the dramatic changes in hydroclimate observed in the geologic record from the Great Basin (and other regions) have motivated substantial work on the response of regional climate to external climate forcings, such as glacial-deglaciation transitions (e.g., Broecker, 2010; Hudson et al., 2017; Ibarra et al., 2014; McGee et al., 2018; Reheis, 1999).

There are multiple hypotheses surrounding the timing and importance of various mechanisms driving changes in LGM and deglacial hydroclimate in the western United States. At present, the northwestern Great Basin’s greatest seasonal precipitation occurs during the winter months. Hence, one group of hypotheses has centered around the LGM and deglacial response of the cool-season mid-latitude jet stream and storm tracks to changing climate forcing (Hostetler and Benson, 1990; McGee et al., 2018; Munroe and Laabs, 2013b); a related hypothesis highlights changes in atmospheric rivers and concomitant changes in atmospheric moisture convergence during the LGM (Lora, 2018; Lora et al., 2017). Field and modeling studies have hypothesized that both the strength and position of the jet stream could be important: it is thought that the Laurentide ice sheet deflected the jet stream south during the LGM, shifting the average storm track and resulting in a tendency for lake highstands to occur along a southeast-northwest trend, modulated by ice sheet retreat during the deglacial (Lyle et al., 2012; McGee et al., 2018; Oster et al., 2015). Finally, there is a series of studies that have explored the correspondence between North Atlantic stadials and the growth and decay of deglacial lakes, with several studies suggesting that the location and intensity of the jet stream acted to modulate local deglacial precipitation rates (Hostetler et al., 1994; McGee et al., 2018; Munroe and Laabs, 2013b).

Previous compilations of lake hydrographs for pluvial lakes in the Great Basin suggest that deglacial lake highstands did not occur simultaneously; instead, changes in regional hydrology may have occurred earlier in the southeast, progressing to the northwest through time (Hudson et al., 2017; Ibarra et al., 2014; Lyle et al., 2012; McGee et al., 2018; Santi et al., 2019). A closer look indicates that this interpretation may be an oversimplification of the process; one study found that a strengthened jet stream during the early deglaciation resulted in higher precipitation along most of the west coast (Lora et al., 2016), while others suggested that alternative mechanisms increased effective moisture, such as intensified summer-season moisture during the deglaciation (Lyle et al., 2012) and changes in evaporation rates (Ibarra et al., 2014).

Additionally, there is proxy-derived evidence of depressed evaporation rates occurring in Lakes Bonneville and Lahontan (Fig. 1 inset) (Kaufman, 2003; Mering, 2015; Mifflin and Wheat, 1979) and in the Great Basin as a whole during the LGM through the early deglacial period (Mifflin and Wheat, 1979; McGee et al., 2018). Climate models simulate LGM increases in both summer and winter effective moisture, or the surplus of precipitation over evaporation (P-E), in the region, but driven by both decreased evaporation and increased precipitation, as controlled by various components of the moisture budget (Lora, 2018). Although there are a number of circulation changes that have been suggested to explain the observed changes in hydroclimate, it has been difficult to robustly test them because of a lack of data constraining evaporation and precipitation rates.

From a mass-balance perspective, lake growth or reduction observed through the dating of shore-zone deposits is achieved primarily via changes in precipitation, evaporation, or a combination of the two (Broecker, 2010; Matsubara and Howard, 2009; McGee et al., 2018). Various proxy evidence from the LGM and early deglacial period (e.g., packrat middens, halite inclusions, glaciers, and pollen) indicate generally cold and wet conditions (Galloway, 1970; Lowenstein et al., 1999; Thompson et al., 1999). While studies have attempted to quantify past evaporation and precipitation rates during this period using lake deposits, invoking either reduced or elevated precipitation rates compared to modern, and reduced evaporation rates compared to modern values (Ibarra et al., 2014; Matsubara and Howard, 2009), there is significant uncertainty associated with these measurements, largely due to a lack of accurate constraints on temperature or on lake-water δ18O. Therefore, values for reconstructed precipitation rates range from less than to more than modern (80%–260% of present values) and evaporation rates span ∼10%–90% of modern values, with cooler conditions by 2–15 °C (Antevs, 1952; Galloway, 1970; Kaufman, 2003; Lowenstein et al., 1999; Matsubara and Howard, 2009; Smith and Street-Perrott, 1983; Thompson et al., 1999).

In this work, we use clumped isotopes, a thermodynamically based tool for estimating carbonate mineral precipitation temperatures (Bernasconi et al., 2018; Eiler, 2007; Ghosh et al., 2007; Schauble et al., 2006), to constrain past temperature and water-isotope changes during both the LGM and deglaciation for Lake Surprise, California, a small pluvial lake located in a hydrologic transition zone between the Great Basin and the Pacific Northwest (Egger et al., 2018; Oster et al., 2015; Wise, 2010). Information on lake-level fluctuations and sediment geochemistry–derived data are combined with different sets of assumptions within a hydrological modeling framework to estimate precipitation and basinwide evaporation rates, producing estimates that are thermodynamically based. We use these results to test hypotheses about the timing and magnitude of hydroclimate changes, and to serve as a point of comparison to published pollen-derived estimates of past hydroclimate during the LGM (Bartlein et al., 2011).

Geologic and Climatic Setting of Surprise Valley, California

Lake Surprise was located in the northwestern Great Basin, along the border of Nevada and California, contained within the modern Surprise Valley (Fig. 1). At its greatest extent, Lake Surprise covered 1366 km2, or ∼36% of its watershed (Ibarra et al., 2014; Reheis, 1999). In contrast to these past wet conditions, modern climate is much drier; average precipitation and lake evaporation rates in Surprise Valley and in southern Oregon (at Goose Lake) are 566 ± 165 mm/yr and 1066 ± 80 mm/yr, respectively (Ibarra et al., 2014; Phillips and Van Denburgh, 1971). As a result, modern Surprise Valley has the same arid climate that characterizes much of the Great Basin; however, due to its proximity to the Pacific Northwest, it receives more precipitation than the southern Great Basin.

Previous Work on Past Hydroclimate in Surprise Valley

Published studies of Lake Surprise have tracked the evolution of the lake levels through time using shore-zone tufa, inferring an abrupt increase in effective moisture leading to a highstand at 16 ka, followed by a much slower decline in lake levels (Egger et al., 2018; Ibarra et al., 2014; Santi et al., 2019). One study used the timing of lake-level fluctuations and carbonate oxygen-isotope measurements to constrain a mass-balance model for precipitation rates, but did not have constraints on temperature or water δ18O (Ibarra et al., 2014). Regional pollen studies also provide nearby constraints on past precipitation rates; pollen data from a locality 4° north of pluvial Lake Surprise support a local decrease in precipitation rates during the LGM (19–23 ka), while several pollen estimates south of Lake Surprise suggest the opposite hydrological trend (Bartlein et al., 2011).

Carbonate Clumped Isotope Thermometry

Carbonate clumped isotope thermometry is a geochemical method to constrain past temperatures that can be applied to sediments (Bernasconi et al., 2018; Ghosh et al., 2007; Huntington and Lechler, 2015; Kele et al., 2015; Petryshyn et al., 2016; Tripati et al., 2010). It is based on the measurement of the overabundance of “clumped” or doubly substituted bonds in carbonate groups of minerals (13C-18O-16O) above their stochastic distributions, which is temperature dependent (Ghosh et al., 2007; Schauble et al., 2006). Gas-source mass spectrometry of CO2 produced through the digestion of carbonate in orthophosphoric acid is used to determine the abundance of the doubly substituted isotopologue with a mass of 47 atomic mass units (amu), and the overabundance of this isotopologue in a sample (relative to a stochastic value) is denoted by Δ47, defined as:
where R is the ration of the CO2 isotopologue of mass 47, 46, or 45 over mass 44 CO2.
The utility of clumped isotope analysis lies in the thermodynamic preference for clumped bonds to occur at certain temperatures; clumping decreases with increased temperature (T), and this trend scales with forumla (T in Kelvin). For example, an isotope exchange reaction that forms clumped bonds:
takes place at equilibrium within a single phase, with lower temperatures favoring a greater abundance of 13C-18O “clumped” bonds (Schauble et al., 2006). Hence, the temperature of carbonate formation can be determined from the Δ47 parameter, without knowledge of the isotopic composition of the fluid in which a given sample formed. Analysis of modern microbialites, tufas, and other types of lacustrine carbonates indicates that this proxy can be robustly used to reconstruct temperature, with growth temperatures typically indicating formation in the summer or spring through fall (Effler and Johnson, 1987; Hren and Sheldon, 2012; Santi et al., 2019; Versteegh et al., 2010).

The comparison of clumped isotope paleoclimate reconstruction to those derived from other methods is valuable, as the clumped proxy is founded on a fundamentally different set of geochemical and conceptual assumptions, and has different sources of uncertainty. For one, because the isotope exchange reactions that form the basis of clumped isotope analysis take place within a single phase, water δ18O can be constrained from the carbonate δ18O; conversely, in other paleoclimate applications, water δ18O is commonly estimated or assumed in the analysis.

The clumped isotope system is dependent on several assumptions that are unique to this proxy system. Although not discussed in detail here, sources of uncertainty related to the Δ47 parameter include temperature-dependent carbonate acid digestion offsets, sample purification, secondary carbonate formation, and the standardization process used to convert raw Δ47 to final Δ47. Discussions of Δ47 parameter uncertainties are included in other publications (Huntington et al., 2009; John and Bowen, 2016; Spencer and Kim, 2015). In light of the fundamental differences between clumped isotope analysis and other proxy systems, independent paleoclimate estimates derived using Δ47 can provide an invaluable point of comparison to existing paleoclimate reconstructions.



We measured the stable and clumped isotope compositions of 35 carbonate (tufa) samples dated to the Late Pleistocene from Surprise Valley, as well as one modern playa carbonate. These samples were described and ages determined using radiocarbon and/or uranium-series measurements in prior publications, and are also described in Table DR1 in the GSA Data Repository1 (Egger et al., 2018; Ibarra et al., 2014; Marion, 2016; Santi et al., 2019).

Isotope Analyses

An extended description of sample preparation for mass spectrometry and lab standardization methods are included in the Data Repository (see footnote 1). We use gas-source mass spectrometry for clumped isotope (Δ47) analysis, which simultaneously also provides carbonate δ13C and δ18O measurements. For each tufa with adequate sample mass, we include at least four replicates, because Δ47 standard error is minimized by increasing the number of sample replicates (Fernandez et al., 2017).

We report Δ47 in the absolute reference frame (Dennis et al., 2011) and standardize it using both solid carbonate and gas standards (Dennis et al., 2011; Tripati et al., 2015, 2010) and using the Brand parameter set (Daëron et al., 2016). Average standard values for the duration of the data analysis are included in Table DR2 (see footnote 1). We remove outliers for replicates whose Δ47 indicates unrealistic water temperatures (≥100 °C or ≤0 °C). We also remove replicates where the standard error introduced by the replicate to the sample set exceeds a specified threshold defined by the long-term reproducibility of standards (Meckler et al., 2014).

For all radiocarbon results, we use the IntCal13 calibration curve (Stuiver and Polach, 1977) to convert conventional radiocarbon ages to calibrated radiocarbon ages, expressed as thousands of years before present, “ka B.P.” We plot the median calibrated probability and the 2σ uncertainty. Samples that were previously dated only by uranium series (Ibarra et al., 2014) were dated by radiocarbon in subsequent publications (Santi et al., 2019; Santi, 2019), and we report radiocarbon results (2σ) here.

Water Temperature

There are numerous clumped isotope calibrations in the literature, which can differ in terms of carbonate material analyzed, reference frame used, or even whether biotic or synthetic carbonates were chosen to construct the calibration. To convert Δ47 to water temperature (Tw), we use the clumped isotope calibration deemed most appropriate for our analysis, because it was created using modern lacustrine tufa, the same carbonate material as our data set, and was also created using the same reference frame as this analysis (Arnold et al., 2018; Dennis et al., 2011):

We compare the slope and intercept of this calibration to those of other clumped isotope calibrations in Table 1.

In the Data Repository (see footnote 1), we carry through our complete set of calculations (for temperature, precipitation rate, and evaporation rate) using two additional calibration equations:

Equation 4 is derived from an array of modern fluvial and lacustrine carbonate types, including tufa, travertine, and microbialites (Arnold et al., 2018), while Equation 5 is chosen because it uses the same carbonate standard reference frame as our data (Bernasconi et al., 2018). Temperature results from alternate calibrations are included in Figure DR1.

Using Equations 35 on each sample Δ47 results in a spread of ∼2–5 °C in absolute temperature for each sample. This range is an overestimate of the uncertainty associated with clumped isotope reconstruction, because not all calibrations are equally suitable for application to our tufa samples (Arnold et al., 2018; Kele et al., 2015). Furthermore, relative changes in temperature through time are less affected by the calibration choice than absolute temperatures (Fig. DR1 [see footnote 1]). Hence, in order to provide the most robust reconstruction of past climate, we employ a calibration equation derived from the same carbonate material as our samples, and base our interpretation on temperature anomalies, rather than solely on absolute water temperatures.

Mean Annual Air Temperature (MAAT)

Hren and Sheldon (2012) created a suite of transfer functions that translate water temperature to mean annual air temperature (MAAT). Each of these transfer functions (April–June, June–August, April–October) assumes a different dominant seasonality of carbonate precipitation; that is, these equations assume that carbonate growth is seasonally biased toward warmer months. Existing literature suggests that not all transfer functions are equally suitable for a given sample set. For example, previous studies of mid–high northern latitude biotic carbonate (Versteegh et al., 2010; Yan et al., 2009) show that most biotic shell growth occurs beginning in the late spring and continuing into the summer. Hren and Sheldon (2012) recommended use of the April–June or April–October transfer function for application to mid–high latitude biotic carbonates, and we choose the April–October transfer function:
because this range has been shown to account for >95% of shell growth (Versteegh et al., 2010). In addition, a recent analysis involving modern fluvial and lacustrine carbonate finds that modern MAAT and full growth season air temperatures (April–October) are better correlated than MAAT and the warmest season of the year (e.g., June–August) (Arnold et al., 2018).

Thus, in choosing Equation 6, we assume that our clumped isotope–derived water temperatures are most representative of April–October lake water temperature. As a point of reference, we note that the average MAAT calculated using the April–June transfer function is 3.5 °C warmer, on average, than when calculated using the April–October transfer function. However, independent of which transfer equation is used; we calculate LGM and deglacial MAATs that are (on average) colder than modern. Comparison of MAAT derived from the April–June and the April–October transfer functions is reported in Figure DR4 (see footnote 1).

When reporting uncertainties associated with MAAT, we add standard error from Tw and error associated with Equation 6 (standard error = 1.89) in quadrature.

Precipitation and Evaporation Modeling

To model the hydroclimate drivers of Lake Surprise, we combine the clumped isotope–constrained precipitation and evaporation modeling approach previously used on Lake Bonneville (Mering, 2015) with an isotope mass-balance model previously used for Lake Surprise (Ibarra et al., 2014), which was modified from earlier work (Jones et al., 2007; Jones and Imbers, 2010). Similar isotope-based mass-balance approaches have been applied to both modern transient and Pleistocene steady-state calculations for other mid-latitude lake systems in the western United States and Europe (Ibarra et al., 2014; Jones et al., 2007), but fundamentally have lacked a thermodynamic constraint on temperature or robust estimates of changes in lake water δ18O.

Lake Evaporation Rate

In this work, we use an equation for lake-based evaporation (EL) (Linacre, 1993) that relies on inputs of latitude (Lat), air temperature (T), dew-point temperature (Td), wind speed (u), and elevation (z), and has been used for previous paleoclimate reconstructions (Ibarra et al., 2014; Jones et al., 2007; Mering, 2015):

For our primary calculations, we assume that u and z have remained constant through time, and that Td is offset a constant amount from temperature, which is reasonable, assuming small changes in relative humidity (RH) (Ibarra et al., 2014; Jones et al., 2007; Linacre, 1993). We assume T is equal to MAAT, but this assumption may bias EL to high values, were the lake frozen over (thus inhibiting evaporation) for a significant amount of each year. In our sensitivity analysis, we explore the effects of allowing u and RH to be altered within a reasonable range (Fig. DR2 [seefootnote 1]).


Currently, LGM climate model boundary conditions do not include pluvial lakes in the western United States (Braconnot et al., 2012; Kagayama et al., 2013), though lakes in the southwest have been implemented in recent Pliocene simulations (Pound et al., 2014). For the most direct comparison between our calculations and the climate model–derived land surface evapotranspiration rate (the evspsbl climate model output variable), we calculate past evapotranspiration (ET) (scaled over land) as a residual of the basin-scale water balance: a function of precipitation (P), the runoff coefficient (krun), and the tributary area (At), as defined in Equation 8:

Modern ET is 528 mm/yr at Lake Surprise (Ibarra et al., 2014). Both P and krun are calculated in our model for precipitation (see Equations 1011 below). We acknowledge that given the lack of pluvial lakes in the boundary conditions of the Paleoclimate Modeling Intercomparison Project, phase 3 (PMIP3), climate models, it is possible that our calculations (including free lake evaporation rates) may be overestimates compared to climate model–derived values. Despite this limitation, PMIP3 provides a relevant comparison to our proxy measurements, as the ensemble allows for a direct comparison of proxy data with an array of model simulations for the same time period as our samples- the LGM.

Weighted Evaporation Rate

Lake evaporation rates are typically higher than basin-scale evaporation rates, because water availability limits the amount of evaporation that can occur from unsaturated surfaces. We derive an equation for weighted evaporation (Ew) that scales free evaporation from the lake surface (EL; Equation 7) and our derived ET from the land surface (Equation 8) by lake area (AL) and tributary area (At), respectively:

EL and Ew for each sample are included in Table DR3 (see footnote 1).

Estimating Lake Area and Basin Hypsometry

The pluvial hydrologic index (HI), forumla, is a physical basin parameter that describes the ratio of lake surface area (AL) to watershed or tributary area (At), and is a primary input in our calculation of precipitation rate (see Equation 10 below). Historically, it has been used as a means to determine the partitioning of rainfall into runoff and evaporation and otherwise approximate past hydroclimate, assuming minimal change in drainage area and a basin’s hypsometry (Ibarra et al., 2018, 2014; Mifflin and Wheat, 1979; Reheis, 1999). We calculate the HI corresponding to each sample elevation as a function of elevation using a hypsometric curve derived from the HydroSHEDS and HydroBASINS data sets (Lehner et al., 2008; Lehner and Grill, 2013).

Lake Precipitation Rate

Beginning with the time-varying water balance and δ18O isotope mass-balance equations for an inward-draining lake and applying the product rule, we derive a function for calculating precipitation rate, modified from equations and derivations in Ibarra et al. (2014), Jones and Imbers (2010), and Steinman and Abbott (2013). The Data Repository (see footnote 1) includes a complete derivation of Equation 10:
This equation includes physical basin parameters (HI) modified by the isotopic composition of lake water (δ18OL), input water (δ18Ow), and evaporating water vapor (δ18Oe). In previous work, a value of krun was assumed (Ibarra et al., 2014); however, modern hydrologic observations suggest a nonlinear response of runoff to changes in precipitation (Broecker, 2010; Greve et al., 2015; Ibarra et al., 2018). To better constrain the relationship between runoff and precipitation, we use the single-parameter formulation for the Budyko curve calibrated for the coterminous United States (Greve et al., 2015):
where Ep is potential evapotranspiration (where we assume that EL = Ep) and ω is the adjustable calibrated Budyko landscape parameter. The use of this Budyko framework in terminal basin hydrologic modeling has been demonstrated in spatially explicit hydrologic modeling, and in similar regional modeling (Ibarra et al., 2018) for Plio-Pleistocene watersheds of the Great Basin, justifying the incorporation of ω into this simplified isotope mass-balance framework.

Given knowledge of EL, basin hypsometry, δ18OL calculated from Δ47, calculation or measurement of δ18Ow and δ18Oe values, and assumptions of ω (variables described below), Equations 10 and 11 can be solved simultaneously for P and krun. Because of the nonlinear nature of both equations, we use a root-finding procedure to solve for the unknowns. This is carried out using the multiroot function in the R package “rootSolve” (Mazzia et al., 2014), which uses a numerical Newton-Raphson method to find the roots of the two equations. Errors are propagated through random draws in the Monte Carlo routine (n = 2500 calculations per sample) by bootstrapping RH, Td, and u, and assuming normal distributions for all input variable values (mean and standard deviation) except for ω, which has a skewed gamma distribution as calibrated for the continental United States (Greve et al., 2015).

Prior to implementing the simultaneous solution to Equations 10 and 11, several model variables need to be determined to populate the equations. We estimate EL using Equation 7. Basin hypsometric curves provide constraints on AL and At. δ18OL is calculated from carbonate δ18O using from Δ47-derived temperature and the temperature-dependent equilibrium fractionation factor (Kim et al., 2007). Meteoric water inputs (δ18Ow) into the lake are constrained from the modern average (–14.57‰ ± 0.6‰, 1σ) (Bowen, 2017), where we assume that source-water effects and temperature effects roughly cancel (Ibarra et al., 2014; Jones et al., 2007). Modern meteoric values for δ18Ow are used because the abundance of measurements allows for a rigorous constraint on this value; furthermore, the modern value is likely similar to that of the LGM, given that many major constraints on water δ18O (e.g., sample elevation, distance from shore, etc.) have not changed appreciably during the Late Pleistocene. Additionally, others have provided constraints on average and standard deviation of modern δ18O, thus constraining the likely range in δ18O for past lake water inputs (Ibarra et al., 2014).

Finally, to calculate δ18Oe, we make the following assumptions:

  • (1) We model evaporating vapor δ18Oe based on the Craig and Gordon (1965) evaporation model (as simplified by Gat, 1996). The kinetic fractionation factor, ε (in per mil), is a simple function of RH: 1000 ln(αkin) ≈ ε = 14.2 × (1 – RH/100), where αkin is the kinetic isotopic fractionation factor based on the ration of Rvapor over Rwater, and R is 18O/16O for vapor and water (Gonfiantini, 1986).

  • (2) The atmospheric vapor above the basin is in equilibrium with the incoming rainwater, which is calculated using the temperature-dependent equilibrium fractionation factor equation (Majoube, 1971).

This approach differs from that of previous studies (Ibarra et al., 2014; Jones et al., 2007), which assumed a kinetic fractionation of αkin = 0.994 for u ≤ 6.8 m/s. In similar work for closed-basin lake modeling (Ibarra et al., 2014; Ibarra and Chamberlain, 2015), the kinetic fractionation factor using the above equation (Gonfiantini, 1986) was found to better approximate the range of possible values (given likely variations in RH), and has been used elsewhere (Horita et al., 2008; Sheldon and Tabor, 2009). Precipitation rates and krun calculated for each sample are included in Table DR3 (see footnote 1).

Potential Sensitivity of Results to Modeling Assumptions

Our isotope mass-balance model relies on several assumptions regarding input parameters. We assume that RH, δ18Ow, u, and ω at the LGM were identical to modern values (Broecker, 2010; Ibarra et al., 2014). In the Data Repository (see footnote 1), we include a model sensitivity analysis (for P and EL) to show the effect of varying each of these four parameters within their 1σ range (Figs. DR2–DR3). Otherwise, the mean modern values of each variable are used: RH = 0.58, u = 1.9 m/s, δ18Ow = –14.57‰, and ω = 2.6.

Statistical Analysis

To understand statistical trends in the data, we use Mann-Kendall trend tests to identify whether our data are monotonically increasing or decreasing during the period of interest. This test is commonly used to interpret environmental, climatological, and hydrological data, and does not require the trend to be linear. Mann-Kendall trend tests have also been used for paleoclimate and long-term climate applications (Futter, 2003; Goossens and Berger, 1987; Tingstad et al., 2014). Results of each trend test are included in Table DR4 (see footnote 1).


Shoreline Geochronology and Carbonate δ18O and δ13C Ratios

A synthesis of Lake Surprise elevations (Fig. 2A) shows a rapid rise in lake levels occurring in <500 yr, culminating in a highstand at 16 kyr B.P. (“16 ka”), as determined using calibrated radiocarbon ages from shore-zone tufas (Futter, 2003; Goossens and Berger, 1987; Tingstad et al., 2014). This rapid rise in lake levels (from 1420 to 1470 m during the LGM, 23–19 ka, to 1545 m at the lake highstand) is seen in several other Late Pleistocene pluvial lakes (labeled on Fig. 1), including Lake Franklin and Lake Lahontan (Munroe and Laabs, 2013b; Santi et al., 2019); however, at Lake Surprise, this rapid transgression is followed by a slow decline over the next ∼5 k.y. There is significant variation in carbonate and water δ18O, occurring on a ∼1000 yr time scale, with quasi-periodic behavior during the LGM (Figs. 2B and 3B, respectively). There are local minima in carbonate and water δ18O at ca. 15–16 ka, coincident with the lake highstand at 16 ka. All reconstructed lake water δ18O values suggest significant evaporative enrichment relative to the isotopic composition of precipitation, which for modern day is –14.57‰ ± 0.6‰ (Bowen, 2017). Between a limited age range, 17.5–14.5 ka, water δ18O increases from its pre-highstand minimum, as confirmed by a Mann-Kendall trend test (n = 18). Conversely, there is little evidence for temporal variability in carbonate δ13C throughout our record (Fig. 2C). A Mann-Kendall trend test (n = 23) indicates that δ13C was stable between 21 and 15 ka (Table DR4 [see footnote 1]).

Clumped Isotope Constraints on Past Hydroclimates

Proxy reconstruction and data analysis requires that tufa samples maintain the same climatic values they contained when they were originally precipitated from the lake water and have remained within a closed system since deposition. Previous reconstructions of lake shorelines indicate that Lake Surprise remained a closed basin, even at its maximum extent (Ibarra et al., 2014; Reheis, 1999). Furthermore, strong positive covariance between carbonate δ13C and carbonate δ18O is a typical feature of closed watersheds (Li and Ku, 1997; Talbot, 1990), and is observed in our data (Fig. 2D). Finally, three scanning electron microscope (SEM) images show that tufa from Surprise Valley is composed of dense calcite, and thus represents primary mineralogies (Fig. DR5 [see footnote 1]). Hence, topographical, geochemical, and crystallographic evidence all indicate that primary climate signals are preserved in our tufa samples.

As presented below, in the lake hydrograph we find no evidence for large precipitation increases leading to the near doubling of lake level recorded by geochronology of the shore-zone tufa (Egger et al., 2018; Ibarra et al., 2014; Santi et al., 2019). Given the steep hypsometric geometry of Surprise Valley, this requires that the lake highstand (ca. 16 ka) at Lake Surprise was due to transient increases in moisture balance (increased precipitation and/or decreased evaporation) via atmospheric teleconnections, rather than long-term steady-state changes (McGee et al., 2018; Munroe and Laabs, 2013a). As such, we feel it is appropriate to make steady-state assumptions in our subsequent isotope mass-balance model (Horton et al., 2016; Talbot, 1990).

Our LGM (19–23 ka) samples indicate an average water temperature of 13.3 ± 0.7 °C (mean and standard error reported throughout text), which corresponds to an average MAAT of 4.5 ± 2.0 °C (Hren and Sheldon, 2012). The offset between modern MAAT from Cedarville, California, near Surprise Valley (9.2 ± 1 °C) and our LGM averaged MAAT (23–19 ka) indicates 4.7 ± 2.2 °C of air warming since the LGM.

To a first order, reconstructed water temperatures appear relatively constant throughout the LGM, with a slight decrease during the early to mid-deglacial period, and a local minimum at ca. 16 ka, roughly coincident with the lake highstand (Fig. 3A). A Mann-Kendall trend test supports a decreasing trend in water temperature between 21 ka to 15 ka (Table DR4 [see footnote 1]). The majority of samples (∼90%) indicate that LGM and deglacial water temperatures were lower than modern water temperature, as estimated by translating modern MAAT at Cedarville to water temperature (Tw = 16.8 ± 2.1 °C, by inverting Equation 6).

In Figures 3C–3D, we report reconstructed precipitation and (weighted and lake) evaporation rates. We show that precipitation rates were slightly above their modern values during the LGM, and stabilized during the deglacial period to near-modern values. A Mann-Kendall trend test indicates that precipitation rates were likely decreasing during the late LGM through the mid-deglacial period (Table DR4 [see footnote 1]). Lake evaporation and basin-scale weighted evaporation rates decreased throughout the LGM and stabilized during the deglacial period (as confirmed using Mann-Kendall trend test; Table DR4), with lake evaporation rates falling below the modern lake evaporation rate near Lake Surprise by 18–20 ka.

Pollen-Derived Estimates of Precipitation

Compilations of proxy data provide an invaluable means to quantify past climate, because each proxy is likely sensitive to different components of the water balance. For example, pollen data are thought to be sensitive to changes in available energy during growing seasons, while lake level fluctuations likely reflect changes in effective moisture, or P-E (Liu et al., 2018). Pollen data have previously been used to provide robust quantitative paleoclimate estimates at both regional and global scales (Bartlein et al., 2011), and have been compared and evaluated against other proxy estimates and results from steady-state model simulations (Liu et al., 2018; Lora et al., 2017; Matsubara and Howard, 2009).

For our analysis, we compare precipitation anomalies derived from clumped isotope analysis to pollen and plant macrofossil precipitation anomalies within the Great Basin (Fig. 4). Bartlein et al. (2011) compiled pollen and plant macrofossils samples that were dated between 19 and 23 ka, the same time frame as this analysis. A subset of the Bartlein samples are also taken from the Great Basin, thus making an ideal point of comparison to Lake Surprise tufa. They calculated precipitation anomalies (LGM minus modern) ranging from –1235 to 721 mm/yr, with a mean of –138 ± 545 mm/yr when averaged across the Great Basin (Bartlein et al., 2011). Trends in pollen precipitation reconstruction becomes even more coherent when comparing Great Basin pollen samples on the north and south sides of the Great Basin. These subpopulations of pollen data lie on opposite sides of a climatological transition zone (“dipole”) indicated in the literature (Egger et al., 2018; Hudson et al., 2019; Oster et al., 2015; Wise, 2010). One of these subpopulations encompasses northern Oregon and southern Washington, with an average precipitation anomaly of –792 ± 540 mm/yr, indicating drier-than-modern LGM conditions in the northwestern Great Basin. The other subpopulation of pollen samples is located in the southeastern Great Basin, and results in an average precipitation anomaly of 141 ± 280 mm/yr, indicative of wetter-than-modern conditions during the LGM. Hence, pollen anomalies demonstrate a spatial pattern of precipitation in the Great Basin, and indicate that hydroclimate varied not just in time, but also in space.

Evaluation of Climate Model Simulations of Hydroclimate Change

We also compare our results to model simulations of past climatic states, including nine steady-state models from the Paleoclimate Modeling Intercomparison Project, phase 3 (PMIP3) (see Figs. 46 for model names), and one transient climate model, Transient Climate Evolution (TraCE) (He, 2011; Liu et al., 2009). For TraCE, we use LGM climatological data averaged between 19 and 20 ka, and modern output averaged between 1984 and 1990 CE (Lora et al., 2016). For PMIP3, we use equilibrium simulations for the LGM (21 ka) and the pre-industrial era (1850 CE). Previous work comparing PMIP precipitation simulations for the LGM indicate a general lack of agreement between models (Lora, 2018; Lora et al., 2017; Oster et al., 2015), and thus the comparison of model results to proxy data offers an opportunity to evaluate model skill.

As a method of visually assessing individual climate model skill, we compare MAAT anomalies calculated using LGM age samples (19–23 ka) and modern weather-station data to climate-model temperature anomalies from LGM and modern simulations (Fig. 5). From clumped isotope analysis, we estimate a 4.7 ± 2.2 °C anomaly between the LGM (4.5 °C) and modern (9.2 °C), implying colder climates during the LGM.

In Figures 4 and 6, we compare precipitation and ET anomalies calculated using LGM samples (19–23 ka) and modern weather-station data to climate-model anomalies, with proxy-derived precipitation and ET anomalies of 22 ± 195 mm/yr and –36.5 ± 183 mm/yr, respectively. We note that TraCE does not directly include output for evaporation. While TraCE does include latent heat flux from the surface, which can be converted to an approximate evaporation by using the latent heat of evaporation, we do not make those conversions here.


Shoreline Geochronology and Carbonate δ18O and δ13C Ratios

The lake hydrograph for Lake Surprise (Fig. 2A) indicates there was a rapid increase in lake levels at 16 ka, suggesting a large and rapid change in effective moisture that is also observed in neighboring pluvial lakes (Egger et al., 2018; Hudson et al., 2019; Ibarra et al., 2014; Munroe and Laabs, 2013a; Santi et al., 2019).

In contrast to its rapid rise, lake regression was more prolonged, which is a notable point of contrast to other Late Pleistocene lakes, including Lake Franklin, Lake Chewaucan, and Lake Lahontan, which surround Lake Surprise during the LGM, as indicated on Figure 1 (Munroe and Laabs, 2013b; Hudson et al., 2019; Santi et al., 2019). We suggest that this gradual decrease in lake levels could be due to the significant depth of Lake Surprise at its maximum extent (∼180 m, versus ∼90 m for nearby Lake Chewaucan; Egger et al., 2018), its high hydrologic index compared to more southerly lakes Franklin, Lahontan, and Bonneville (Santi et al., 2019), or the relative lack of western-boundary orographic barriers compared to other lake basins (e.g., the southern Cascades to the west of Lake Chewaucan; Egger et al., 2018).

Clumped Isotope Constraints on Past Hydroclimates

Using modern climate averages from Cedarville, California, and LGM samples from Lake Surprise (23–19 ka), we estimate 4.7 ± 2.2 °C of air warming since the LGM at Lake Surprise. As a point of comparison, Δ47-derived water temperatures from shore-zone tufa at nearby Lake Chewaucan for the modern day (13 ± 2 °C) and late LGM (ca. 18 ka; 6.2 ± 2 °C) indicate 6.8 ± 2.8 °C of lake water warming, and 10.0 ± 3.3 °C of air warming (Hren and Sheldon, 2012; Hudson et al., 2017).

Our LGM anomalies are derived from 17 individual samples dated between 23 and 19 ka spanning ∼30 m of shoreline elevation (Table DR1 [see footnote 1]), and include samples from the majority of the shoreline localities sampled by Egger et al. (2018) and Ibarra et al. (2014), and thus correspond to a range in clumped and stable isotope values. The LGM samples have a Δ47 standard deviation of 0.018‰ (0.004‰ standard error), which is similar to the predicted variability in Δ47 associated with analytical error (Fernandez et al., 2017; Huntington et al., 2009) and corresponds to a temperature variability of ±4 °C. With respect to stable isotope values, the sample average (μ) and standard deviation (σ) for δ13C (μ = 3.5‰; σ = 0.25‰) and δ18O (μ = –3.4‰; σ = 0.29‰) show limited variability, and are furthermore in good agreement with values from LGM tufa measured at Lake Chewaucan in southern Oregon (Hudson et al., 2017). Finally, the slight rise in lake level recorded by the sample elevations is likely due to small transient increases in the water balance over the LGM, though we do not see corresponding changes in δ18O. Because Surprise Valley has a steep hypsometry at these elevations (i.e., a minimal lake area change for substantial changes in elevation), this suggests a similar steady-state water balance over the entire LGM period.

Many paleoclimate studies have produced or compiled estimates of past climate at larger spatial scales than this work, including in the Wasatch Range, located on the eastern edge of the Great Basin (Quirk et al., 2018), the Mohave Desert in the southern Great Basin (Quade et al., 2003), and in the Great Basin as a whole (Mifflin and Wheat, 1979). These studies have produced estimates of temperature depression ranging from 2 to 11 °C. Several of these suggest greater temperature depressions than implied by clumped isotope analysis, including ones based on pollen (10–11 °C; Galloway, 1970), hydrologic mass-balance modeling (10 °C; McGee et al., 2018), glacial modeling (8–10 °C; Quirk et al., 2018), and packrat midden plant assemblages (8 °C; Thompson et al., 1999). Other Great Basin–scale analyses are more in line with our results; two mass-balance models suggest 2–3 °C of air cooling in the Great Basin during the LGM and deglacial period (Antevs, 1952; Mifflin and Wheat, 1979). Comparison to other LGM and deglacial paleoclimate reconstructions must also acknowledge the use of overlapping but unidentical time frames for analysis, including 23–17 ka (Galloway, 1970), 26–16 ka (McGee et al., 2018; Quirk et al., 2018), and 20.5–18 ka (Thompson et al., 1999). Through this analysis, we define the LGM as the period from 23 to 19 ka, as this encompasses the PMIP LGM time slice, and modeling studies suggest that 19 ka is when the Laurentide Ice Sheet began its retreat across the North American continent and greenhouse gas concentrations began to rise (Lora, 2018; Lora and Ibarra, 2019; Marcott et al., 2014; Peltier, 2004). Despite variation in the time periods that are represented in previous work, we note that our paleoclimate reconstructions fall within the range of proxy estimates.

Further, as discussed in Lora and Ibarra (2019), some of the inconsistencies between paleoclimate data and the existing deglacial climate modeling experiment (TRaCE) may be due to the timing of ice sheet elevation and extent, and the timing and implementation of meltwater pulses (Deschamps et al., 2012; Ivanovic et al., 2018; Menounos et al., 2017). As new ice sheet reconstructions are being implemented (Peltier et al., 2015; Tarasov et al., 2012) into the next generation of LGM (PMIP4) and deglacial climate-modeling experiments (Ivanovic et al., 2016), we expect these inconsistencies to be resolved, and for future models to simulate transient increases in precipitation (e.g., McGee et al., 2018), which may have driven the brief lake highstands established during the deglaciation, as observed here and further inland.

A more direct comparison to this work is the modeling study for Lake Surprise by Ibarra et al. (2014). They derived past precipitation and evaporation rates based on analysis of shore-zone tufas (some of which were analyzed in this study) from Lake Surprise, assuming 7 °C of air warming since the LGM period and 5 °C of warming between the deglacial period and modern day. Their assumption was based on combining the stable isotope, elevation, geochronologic, and hypsometric data from Lake Surprise with nearby pollen-assemblage data and temperature depression estimates from Little Lake, Oregon, to the northwest (Worona and Whitlock, 1995). Significantly, our estimate is also in line with a global synthesis of pollen data, which indicates 4.4 ± 0.9 °C of air cooling during the LGM at sites in the vicinity of Lake Surprise (Bartlein et al., 2011).

Variation in water δ18O has been traditionally used to imply variation in paleotemperature at the site of water condensation (Dansgaard, 1964), but has more recently been used to also glean information on moisture source, seasonality of precipitation, and rainout history (Charles et al., 1994; McKenzie and Hollander, 1993; Plummer, 1993). Qualitatively, our samples indicate fairly constant water δ18O with a local minimum between the late LGM (ca. 19 ka) and lake highstand (ca. 16 ka) (Fig. 3B). This depletion could be due to changes in the dominant lake moisture source; for example, the δ18O of water associated with the North Pacific storm track is depleted in heavy isotopes relative to the δ18O of water associated with atmospheric rivers (Welker, 2012). We note that summer precipitation in the western United States is enriched in heavy isotopes relative to winter precipitation, so this pattern could reflect an increase in winter precipitation occurring near the time of the lake highstand compared to the late LGM (Welker, 2012). A decrease in water δ18O could alternatively be explained by decreasing condensation temperature, with lake temperature fundamentally related to the temperature of the overlying air (Dansgaard, 1964). Finally, this variability could be due to changing levels of lake evaporation through time.

To a first order, evaporation rates at Lake Surprise decrease throughout the LGM and deglaciation, with lake evaporation rates stabilizing below the modern lake evaporation rate of 1066 mm/yr by ca. 18–20 ka (Fig. 3D). Weighted evaporation rates are consistently lower than lake evaporation rates. This is expected because weighted evaporation is a weighted average of lake evaporation rate and evapotranspiration rate, scaled by lake and tributary area respectively, where evapotranspiration rates are lower than free-water lake evaporation.

There are several potential sources of uncertainty in any equation used to estimate evaporation rates. Empirical equations for evaporation are typically reliant on three main categories of controlling parameters: water supply, energy for evaporation, and water vapor transport. The first of these, water supply, can limit evaporation rates in water-scarce or arid regions, or in locations where the evaporating surfaces are frozen over. If Lake Surprise was frozen for significant portions of each year, actual evaporation rates would have been lower than our results suggest. However, because our calculated lake evaporation rates are generally lower than modern lake evaporation rates during the deglacial period, our conclusions would not be significantly altered if we are overestimating lake evaporation.

Our formulation for lake evaporation instead incorporates measurements constraining the latter two parameters: energy for evaporation, and water vapor transport (Linacre, 1993). There is a variety of methods to constrain the amount of available energy to drive evaporation. We choose Equation 7 (Linacre, 1993) because it uses parameters that can either be recorded for each carbonate sample upon collection (elevation and latitude) or constrained through geochemical analysis (temperature). Water vapor transport is also included in many models for evaporation, as it provides a mechanism to remove saturated air from above the evaporating surface. Many evaporation models incorporate measurements of wind speed (u), relative humidity (RH), or vapor pressure deficit into evaporation estimates (Allen et al., 1998; Hargreaves and Samani, 1986; Jensen and Haise, 1963; Penman, 1948; Priestley and Taylor, 1972). Our lake evaporation model incorporates estimates of u, as this is commonly measured at weather stations across the United States. As past wind speed is not a straightforward variable to constrain, we use the modern average value at Surprise Valley, but include this parameter in our sensitivity analysis in the Data Repository (Fig. DR2 [see footnote 1]). Although modest increases in LGM wind speed have been modeled over the ocean off the west coast of the United States (e.g., Lora et al., 2016), these changes have not been found over land near Lake Surprise, at least not on a resolvable spatial scale.

Water vapor transport is also described in Equation 7 using estimates of dew-point temperature, which is a function of RH. We assume modern values of RH in our calculations for past evaporation rates; however, RH is also included as a parameter in our sensitivity analysis. We note that elevated RH would cause dew-point temperature to fall closer to air temperature, resulting in lower calculated past evaporation rates. Hence, even if LGM relative humidity were higher than modern values, our conclusion (lowered evaporation rates) would remain the same.

In Figure 3C, we show estimates of past precipitation rates using our isotope mass-balance model. During the LGM, we find precipitation rates to be elevated (but variable) relative to modern; however, by the time of the lake highstand at 16 ka, calculated precipitation falls to near-modern values. This finding is significant, as it implies that lake growth (and increased effective moisture) was achieved despite LGM precipitation rates that were similar to modern. While the rapid lake growth at 16 ka suggested by our data (106 m in ∼500 yr) does require a significant short-term surplus of precipitation over evaporation, the short time scale over which lake transgression occurred (Fig. 2A) coupled with the comparatively coarse sampling resolution of our data could explain why this positive P-E anomaly is not fully reflected in Figure 3C.

Our sensitivity analysis shows that our model for precipitation rate is most sensitive to changes in RH and u. Again, we note that no significant increases in surface wind speed have been modeled at the location of Lake Surprise (Lora et al., 2016). However, increased RH would result in lowered reconstructed precipitation rates, which remains a source of uncertainty in this analysis.

Comparison with Pollen and Other Proxy Data

Existing Great Basin proxy estimates of precipitation and evaporation rates for the LGM and deglacial have been synthesized in other sources (Lora and Ibarra, 2019; Matsubara and Howard, 2009; Oster et al., 2015). These analyses include climatological averages that represent a considerable temporal range; for example, our work references proxy reconstruction or modeling studies representative of 26.6 ka to 13.5 ka (e.g., Bartlein et al., 2011; Hostetler et al., 1994; Lemons et al., 1996; Quirk et al., 2018). Importantly, our proxy evidence has a relatively high temporal resolution, which could allow for more explicit comparison to other proxy systems in the future.

Our thermodynamically derived LGM-averaged precipitation rate (588 mm/yr) is in line with several Great Basin estimates calculated using different types of proxies (80%–130% of modern) during the LGM and deglacial period (Galloway, 1970; Hostetler et al., 1994; Lemons et al., 1996; Thompson et al., 1999). Figure 4 shows a comparison of clumped isotope–derived precipitation anomalies with an extensive pollen data set across the Great Basin. These pollen precipitation anomalies indicate both wetter and drier LGM conditions, with generally drier conditions in the northwest and wetter conditions in the southeast (Bartlein et al., 2011). We calculate a small positive precipitation anomaly at Lake Surprise, indicating slightly more precipitation during the LGM. Thus, our estimates are consistent with pollen estimates, which suggests a transition from drier-than-modern LGM conditions to wetter-than-modern conditions along a northwest-southeast trend.

As shown in Figure 3D, most EL and Ew are below the modern pan evaporation rate in Surprise Valley (905 mm/yr; Ibarra et al., 2014), as well as the modern EL at Goose Lake near Lake Surprise (1066 mm/yr), falling below the modern average by 18–20 ka (Phillips and Van Denburgh, 1971). We use statistical analysis to show that both EL and Ew decreased beginning in the LGM, continuing into the early to mid-deglacial periods (Table DR4 [see footnote 1]), indicative of a significantly altered climatic state. Our evaporation rates are consistent with paleoclimate indicators from similar time periods, including mass-balance models (–10% relative to modern day; Mifflin and Wheat, 1979), thermal evaporation models (–42%; Hostetler and Benson, 1990), and pollen (–50%; Galloway, 1970), providing evidence that low evaporation was critical to growing and maintaining Lake Surprise during the LGM and deglacial period.

Evaluation of Climate Model Simulations of Hydroclimate Change

Figure 5 shows surface temperature anomalies between the LGM (21 ka) and the modern day, as simulated by PMIP3 and TraCE. We overlay MAAT anomalies calculated using an average from all LGM samples (19–23 ka) and modern climate data. We find reasonable agreement between climate model and proxy estimates (±2 °C) for all models except FGOALS, with all models and proxy data suggesting lowered air temperatures during the LGM.

In Figure 4, we plot precipitation anomalies indicated by climate model simulations and proxy data. Using our LGM samples from Lake Surprise (19–23 ka), we calculate a small positive precipitation anomaly, suggesting only slightly higher precipitation rates during the LGM. This finding is significant, as most PMIP3 climate models show a transition from wetter-than-modern to drier-than-modern climates along a line that is projected through northern california. The exact latitude of this transition zone varies between models; CNRM-CM5 shows this transition along the California-Nevada border at Lake Tahoe, while NCAR-CCSM4 shows this transition occurring along the border at the approximate location of Lake Surprise. With respect to our reconstructed LGM precipitation rates, we infer that Lake Surprise was located near this transition zone (Lora et al., 2016; Oster et al., 2015). Visually, we find reasonable model-data agreement (±100 mm/yr) with all models, as model-data discrepancies are all within error of our clumped isotope–derived estimate.

In Figure 6, we overlay the average ET anomaly indicated by climate model simulations and proxy data. We subtract the LGM-average ET [calculated as P × (1 − krun)] from modern ET at Lake Surprise. We find that all models agree on the sign of the evaporation anomaly (negative), but models MIROC-ESM and MRI-CGCM3 model greater evaporation depression than our samples. Overall, models CNRM-CM5, COSMOS-ASO, and GISS-E2-R (p150) show the best qualitative agreement with proxy reconstructions of T, P, and E using clumped isotopes, indicating moderate elevated precipitation rates and lowered evaporation rates during the LGM.


Constraining the timing and drivers of lake highstands has important implications for understanding the terrestrial and atmospheric processes that transport moisture and impart changes on the basin-scale hydrological cycle. In the Great Basin, paleoclimate studies provide insight into how this region may evolve in the future, as global warming imparts changes on regional air temperatures. In this study, we first use previously compiled lake-level histories to summarize major events in the growth and decay of the Late Pleistocene Lake Surprise. We see evidence of a fast lake transgression at ca. 16 ka (from 1420–1470 m during the LGM to 1545 m at the lake highstand), suggesting a large effective moisture forcing. While the abruptness of this lake transgression is similar to that of nearby Lake Lahontan and Lake Chewaucan, the regression of Lake Surprise occurs more gradually than in these lake basins. We suggest that this slow regression may be due to lake geometry, with both lake depth and hydrologic index playing a role, or due to the smaller orographic barriers to the west of Lake Surprise compared to other lakes, like Lake Chewaucan.

Stable isotope data suggest a minimum in δ18O of precipitation, coincident with the lake highstand at 16 ka; we note that this trend could be explained by changing water sources, changing paleotemperature at the site of condensation, and/or variability in lake evaporation rate.

In addition to lake shoreline data, estimates of past temperature from clumped isotope analysis on lake sediments offer further insight into past hydroclimate. We calculate an average surface water temperature of 13.3 ± 0.7 °C from our LGM age samples, which translates to a MAAT of 4.6 ± 2.0 °C. Using modern air temperature at Cedarville, California, we estimate 4.7 ± 2.2 °C of air warming since the LGM. From Δ47-derived temperatures, we estimate past evaporation rates and precipitation rates. The growth of Lake Surprise implies that the lake hydrologic budget briefly departed from steady state during the deglacial period. Based on our precipitation and evaporation reconstruction, we infer that this departure and subsequent lake growth could have been aided by both decreasing evaporation rates and increasing precipitation following the LGM.

We perform qualitative assessments of model skill at the location of Lake Surprise. We use clumped isotope–derived anomalies of temperature, precipitation, and evapotranspiration on PMIP anomaly plots to visually assess which models are best able to reproduce the hydroclimate anomalies implied by our data and pollen-derived estimates of past climate. Although this first-order assessment only allows for a broad comparison of clumped isotope–derived values to climate models, further compilation of proxy data across a spatial transect will allow for a more quantitative assessment of model skill, which could be used to refine climate models and eventually inform climate policy.

Ultimately, this work sheds light on drivers that supported ancient large-scale lakes in the western United States, and why they disappeared, representing a proof of concept for a method that is broadly applicable to paleoclimate reconstructions and model evaluation using sediments from small closed-basin lakes. The expansion of this approach to multiple lake sites across Western North America, in concert with isotope-enabled simulations, will improve confidence in understanding of regional hydrology, by enabling better determination of past temperatures, evaporation rates, and water vapor sources.


We thank the Tripati Lab Group at University of California, Los Angeles (Team Star Wars) for their support, as well as Anne Egger and Brian Marion for sharing samples dated in Egger et al. (2018) and for sharing their knowledge of Surprise Valley’s Quaternary history, Sarah Lummis for assistance in the field, and Jessica Oster for detailed discussions. We thank Matthew Coble and Conni De Masi for help with SEM imaging. We acknowledge the modeling groups the Program for Climate Model Diagnosis and Intercomparison (PCMDI) and the World Climate Research Program’s Working Group on Coupled Modeling (WGCM) for their roles in making available the WCRP CMIP5 multi-model data set. Finally, we would like to thank Brad Singer and Benjamin Laabs for their careful reading of the manuscript and for prviding suggestions and constructive feedback.

Funding: This project and all University of California, Los Angeles participants were supported by a U.S. National Science Foundation (NSF) CAREER award (NSF EAR-1352212), with mass spectrometry supported by the U.S. Department of Energy through Basic Energy Sciences grant DE-FG02-13ER16402 to AT. LMS and AJA received support from the UCLA Center for Diverse Leadership in Science, and AJA was also supported by a Cota-Robles Fellowship. Kate Maher (Stanford University Earth Systems Science) provided funding for Lake Surprise sample collection by DEI and CW, supported by NSF grant EAR-0921134. DEI was supported by a Heising-Simons Foundation grant to C. Page Chamberlain, the University of California President’s Postdoctoral Fellowship, and the Miller Institute Research Fellowship. JML was supported by University of California Chancellor’s and California Alliance Postdoctoral Fellowships.
Author contributions: LMS led the project and drafted the manuscript. The project was developed by AT with input from JAM and DEI. AJA, CW, DEI, and RL helped collect, process, and analyze samples. Clumped isotope analyses were done by LMS, AJA, JAM, and RL, with the supervision of AT. DEI developed the isotope mass-balance model, and was the driving force in data collection. Hydrologic modeling was done by LMS with input from DEI, AJA, JAM, and AT. JML contributed to the climate modeling component of the work. All authors contributed to editing the manuscript.
Competing interests: There are no acknowledged conflicts of interest.
1GSA Data Repository item 2020190, analytical methods, discussion of the precipitation and evaporation modeling assumptions, supporting figures, and tables, is available athttp://www.geosociety.org/datarepository/2020 or by request to editing@geosociety.org.
Science Editor: Bradley S. Singer
Associate Editor: Benjamin Laabs
Gold Open Access: This paper is published under the terms of the CC-BY license.