Quantifying closed-basin lake temperature and hydrology by inversion of oxygen isotope and trace element paleoclimate records

Lake systems are important paleoclimate archives that preserve ecosystem and hydrologic responses to critical periods in Earth history, such as carbon cycle perturbations and glacial-interglacial cycles. Geochemical measurements of biogenic carbonate (for example, δ18O, δ13C, 87Sr/86Sr, [Li], [U], [Sr], and [Mg]) are indicators of hydrologic variability in lake systems throughout the geologic record. In this study, we present a new closed-basin lake modeling approach, HyBIM (the Hydrologic Balance Inverse Model) that employs a system of total differential equations and uses the measured δ18O, Sr/Ca, and Mg/Ca of biogenic carbonate to determine changes in temperature, runoff, and lake evaporation. Using equally-spaced time steps, these equations are simultaneously solved to constrain the hydrologic parameters of the lake as recorded in biogenic carbonate. We use a Monte Carlo approach to account for uncertainty in the input parameters, such as δ18O temperature relationships, partition coefficient uncertainty, and watershed solute chemistry. For illustrative purposes, we apply the model to two ostracod valve datasets covering different timescales: (1) the Cretaceous Songliao Basin, northeast China, and (2) Holocene Lake Miragoane, Haiti. Modern water measurements of water isotopes and cation concentrations from each location are required as model inputs. We compare our modeling results with author interpretations and geologic observations. The modeling approach presented in this study can be applied to other closed-basin lake records, can be modified for other calcifying species (for example, gastropods or mollusks) or with calibration to inorganic lacustrine carbonate. In addition, this approach holds promise for extension with additional proxy measurements (that is, δD, U/Ca or Li/Ca) and changing source area on tectonic timescales using proxies that reflect changing source lithology (that is, Sr and Pb isotopes). Future incorporation of age model uncertainty in the Monte Carlo approach will also provide utility by quantifying temporal uncertainty on the hydrologic response recorded by lake sediments.


INTRODUCTION 40
Lake systems provide valuable paleoclimate archives as natural integrators of terrestrial 41 processes. Lake deposits from hydrologically closed basins or terminal lakes record important 42 changes in climate, human activity, hydrology, ecology, and tectonics (Antevs, 1948;Eugster 43 and Jones, 1979; Currey, 1990 open versus closed hydrology of a lake system using facies identification has been extensively Horton and others, 2015). We do not incorporate variations in δ 13 C in this initial iteration of 126 HyBIM; however, evaporative processes and lake residence times should also control the 127 systematics of the δ 13 C of carbonates from terminal lake systems (see discussion in Horton and 128 others, 2015). 129 First, we establish the modeling framework using a system of total differential equations. 130 Second, we describe the Monte Carlo routine used to account for parameter and input data 131 uncertainty (for example, Fantle, 2010; Royer and others, 2014). Third, we discuss the input 132 dataset (the time series of measurements) and necessary input parameters. Finally, we describe 133 the quantification of the partial derivatives necessary to constrain the model and necessary input 134 parameters. This model is available as a generic R code built for input of data from any time 135 series of biogenic carbonate from the authors website (http://paleoclimate.stanford.edu). For 136 illustrative purposes, we apply HyBIM to two ostracod valve datasets covering different 137 timescales: (1) the Cretaceous Songliao Basin, northeast China, and (2) Holocene Lake 138 Miragoane, Haiti. All model parameters are summarized in table 1 and a schematic of the model 139 structure is presented in figure A1. 140

Model Framework 141
Here, we present a modeling framework that produces a time series of temperature 142 changes, effects of evaporation, and relative volume estimates for runoff for carbonate-bearing 143 (here ostracods) paleolake records. Chivas and others (1986a) previously suggested a 144 quantitative framework by which a unique solution for paleotemperature and paleosalinity could 145 be calculated using paired Mg/Ca and δ 18 O measurements. Following this observation, we write 146 three total differential equations for lake water that describe the effects of temperature, 147 evaporation, and runoff/precipitation on the oxygen isotopic and trace element concentration of 148 ostracods in a hydrologically closed lake. The three total differential equations (for lake water 149 δ 18 O, [Mg], and [Sr]) are: 150 where δ 18 O is the oxygen isotopic composition of the lake water determined from 154 biogenic carbonate measurements corrected for species-specific vital effects, [Mg] and [Sr] are 155 the lake water concentrations of Mg 2+ and Sr 2+ determined from the Sr/Ca and Mg/Ca ratios 156 measured in the biogenic carbonate, T is mean annual surface air temperature (K) for the lake 157 basin, F evap is the fraction of the lake evaporation, and F input is the fractional lake volume 158 increase due to changes in runoff and precipitation (see table 1). 159 We treat equations 1-3 as a system of total differential equations that can be solved for where matrix A is composed of the partial derivatives, the b vector is the first derivative of the 162 input dataset (left side of eqs 1 to 3) and x is the solution vector of (dT, dF evap , and dF input )). We 163 do so by determining the partial derivatives for a specific lake system, calibrated using modern 164 watershed/lake chemistry data and knowledge of the specific ostracod species. Interpolation to 165 an evenly spaced time series is required, which we accomplish using an Epanechnikov kernel 166 smoother applied to the input data. We assume that within each time step, lake evaporation 167 and/or lake volume change does not change more than 20%, and we approximate all derivatives 168 according to this assumption (that is at 10%). The model output is a time series of the solution vector (dT, dF evap , and dF input ). dF evap and dF input are combined to give the fractional net volume 170 increase for each time-step (see fig. A1). We use a Monte Carlo approach to account for 171 uncertainty in the input parameters and in the input dataset, which we derive from kernel 172 smoothing. 173 174

Monte Carlo routine for uncertainty quantification 175
Given the assumptions that underlie the input parameters and data set for this model it is 176 critical to quantitatively assess the errors associated with the calculated runoff input, temperature 177 changes, and evaporative effects. Thus, we use a Monte Carlo routine to account for two types 178 of uncertainty. The first is uncertainty in the input dataset. Since most paleoclimate records 179 display uneven variance and non-uniform sampling density, it is necessary to account for this by 180 producing datasets for each iteration based on the statistics (time series of the mean and 181 residuals) derived from kernel smoothing of the dataset (see below). To do so, for each Monte 182 Carlo iteration we implement a bootstrapping method whereby we account for input dataset 183 uncertainty by re-sampling residuals. The second type of uncertainty accounted for by our model 184 is the uncertainty in the input parameters. The input parameters are best determined by modern 185 watershed and lake chemistry (when available) and knowledge of the species precipitating the 186 biogenic carbonate (described below). Within each Monte Carlo iteration, we derive a time series 187 of partial derivatives based on the input parameters. Input parameter uncertainty distributions are 188 assumed to be normal distributions, requiring that uncertainty is quantified by as many modern 189 measurements of watershed chemistry as possible. The partial derivatives are substituted into the 190 A matrix at each time step. 191

Input dataset requirements and interpolation 193
There are two requirements for the input dataset that need to be met for this model.  For lakes with no modern analogue (see Cretaceous Songliao Basin example), first-order 212 estimates for τ res can be approximated using scaling relationships between outcrop size and basin 213 geometry (Hendriks and others, 2012), which in combination with facies identification can be 214 used to determine estimates of the relationship between potential accommodation and basin should be used with a time series produced by a single ostracod species or genera. Although, if 239 multiple species are measured (for example, Lister and others, 1991) this could be accounted for 240 using vital effect offsets. 241 In contrast, correcting for temperature effects is less straightforward and could, if done 242 incorrectly, lead to "circular reasoning" since we are attempting to constrain temperature 243 changes in the model. We assume temperature could vary from 10 to 35 °C (uniform 244 distribution). For each Monte Carlo iteration we select a starting temperature and assume that 245 ). 272

Temperature Partial Derivatives 273
To solve for the oxygen isotope partial derivative requires knowing how δ 18 O of 274 meteoric water (precipitation and runoff) varies as a function of temperature. We recognize that 275 there are numerous factors affecting the δ 18 O of precipitation beyond that of temperature, which 276 include seasonality, vapor recycling, moisture source et cetera. However, the commonly held 277 assumption that underpins terrestrial paleoclimate studies is the well-known positive correlation 278 between δ 18 O values of precipitation and temperature for temperate and high latitude areas 279 (Dansgaard, 1964). For the mid-latitude site (Songliao Basin) we use the observed relationship 280 between temperature and precipitation Rozanksi and others (1993) of 0.58 ‰/K, as has been 281 used in numerous studies to calculate the temperature for paleoclimate archives in mid-to high-latitude sites, such as paleosols (for example, Koch and others, 2003), paleolakes (for example, 283

Anderson and others, 2001), and ice cores (see Jouzel and others, 1997 and references therein). 284
For the low latitude site (Lake Miragoane, Haiti) it is unlikely that the correlations of Rozanksi 285 and others (1993) apply as these are for mid-to high-latitudes. Although it is not ideal to use the 286 approach we outline in this paper to low latitudes because temperature is not as well correlated to 287 the δ 18 O of precipitation, we use this site only as an example of how this approach can be used in 288 lakes because this record offers the type of data that are necessary for these calculations. In our model, thus, one of the assumptions inherent in our calculations is that temperature 293 is the primary driver of δ 18 O value of water input to the lake. Note that since the slope between 294 temperature and δ 18 O of precipitation can and does vary geographically, more accurate 295 approaches can the tailored for individual sites by measurements of site-specific correlations 296 between temperature and δ 18 O of precipitation. 297 In addition, to the partial discussed above it is necessary to place the ostracod 298 values as the δ 18 O of lake water for these calculations (b vector of eqs 1 to 3). To do this we need 299 to account for both vital effects and the temperature-dependent isotopic fractionation (see above). 300 In addition, there are minor temperature effects on the Sr/Ca (typically none) and Mg/Ca 301 (see above discussion). Since the inputs (b vector of eqs 1 to 3) are placed in terms of lake 302 composition, after accounting for the Mg/Ca or Sr/Ca of the lake water and selected temperature carbonates, we apply the applicable distribution coefficient (K D ) and take the first derivative. In 312 doing so, we assume that the Sr/Ca and Mg/Ca measured in the biogenic carbonates are faithfully 313 recording changes in the Sr and Mg concentrations of the lake water. Assuming mass 314 conservation, the change in solute M concentration (from the initial concentration to an increased 315 concentration at time 2) during fractional lake evaporation (F evap , fraction remaining, from 0 to 1, 316 where 1 is fully evaporated) is given by the equation (fig. 1): 317 Solving for the partial derivative of equation 5 has the analytical solution: 319 This relationship is non-linear; thus, to provide a linear slope for ∂[M] ∂F evap we assume that for any 321 given time step, the lake does not evaporate more than 20% by volume. By inspection of 322 equation 6, the average evaporation between 0 to 20%, given by F evap = 0.1, yields a slope of parameters are still required, making the generic application of these lake water evaporation 334 models to paleoclimate records challenging. 335 Central to the difficulties encountered by previous researchers is the influence of 336 humidity on the evolution of the isotopic composition of an evaporating water body (Gonfiantini, 337 1986). This is due to changes in the kinetic enrichment factor due to changes in humidity (Gat, 338 1970;Merlivat, 1978;Merlivat and Jouzel, 1979). However, if within a given time step lake 339 evaporation does not exceed ~25% total lake volume and humidity is <90%, the isotopic 340 evolution of an evaporating water body can be approximated by a Rayleigh relationship 341 (Gonfiantini, 1986). Thus, we derive the partial derivative for δ 18 O and fractional lake 342 evaporation as the first derivative of the Rayleigh equation ( fig. 1): 343 where α is the combined fractionation factor of the temperature-dependent equilibrium 345 fractionation factor (Horita and Weslowski, 1996) and humidity-dependent kinetic fractionation 346 factor (Gonfiantini, 1986), F evap is the fraction of the lake evaporation (as above) and δ 18 O initial is 347 the initial isotopic composition of the lake from the measured δ 18 O of the biogenic carbonate 348 corrected for vital effects (see below). While a simplification of evaporative processes, assuming 349 a Rayleigh relationship and a combined fractionation factor in this manner has successfully been 350 applied to both surface soil reservoir and lake system modeling (for example, Chamberlain and 351 others, 2014; Caves and others, 2015). We assume that relative humidity for each Monte Carlo 352 iteration can vary from 50 to 90 % (uniform distribution), and temperature varies from 5 to 35 °C 353 (uniform distribution) as determined for each iteration (as described above). By doing so we 354 calculate the combined fractionation factor, α, using the equations of Horita and Weslowski 355 (1996) and Gonfiantini (1986) (see fig. 1). The humidity and temperature ranges selected here 356 are conservative and broadly apply to both illustrative lake systems used in the examples below. 357 The relevant latitude and geographic setting of the lake system being modeled should inform the 358 humidity and temperature ranges chosen as model input. 359 360 Lake Input Partial Derivatives 361 Determination of the lake input partial derivatives is similar for both trace elements and 362 oxygen isotopes. We rely on binary mixing equations assuming that the isotopic and chemical 363 composition of water entering the lake is sufficiently well characterized and relatively invariant 364 (relative to the lake water). Water in terminal lakes is sufficiently evaporatively enriched relative 365 to the input δ 18 O that in most examples, the input water δ 18 O composition is typically statistically 366 distinguishable from the lake water δ 18 O composition (see for example compilation by Horton and Oze, 2012). This is similar for trace elements in terminal lake systems, such that the input 368 concentration of Sr and Mg in terminal lakes is typically more dilute than the lake water. 369 Given these principles, we define the partial derivative relating increased lake volume, 370 from the addition of source water (F input , which varies from 1 to infinity, where 1 is the original 371 lake volume) and variations in δ 18 O as the first derivative of a binary mixing relationship ( fig. 2): 372 where δ 18 O source is the meteoric water composition and δ 18 O initial is determined as previously 374 This can best be done by measuring the elemental and isotopic composition of modern rivers and 388 precipitation in the lake to be studied or the site of the ancient lake as is done here in the 389

Ibarra and Chamberlain -19
Songliao Basin (table 2). There are obvious complications to this approach since the drainage 390 basin may have changed through time as well as the isotopic composition of precipitation -for 391 example by a change in moisture source. In addition, the uncertainty (standard deviation) of the 392 end-members is required. This can be quantified through repeat measurements of stream 393 water/subsurface runoff chemistry. 394 Second, empirical relationships for temperature effects on the species-specific partition 395 coefficient (for Mg/Ca) are required; such as those fit by Engstrom and Nelson (1991) and 396 DeDeckker and others (1999). As discussed above, species-specific vital effects for δ 18

EXAMPLE APPLICATIONS 406
We apply HyBIM to two lacustrine paleoclimate records covering differing timescales, 407 geochemical conditions, ostracod species, and dataset resolutions (figs. 3 and 4). We chose these 408 two data sets because: (1) they had the elemental and isotopic measurements necessary for using 409 this model and; (2) they represent two very different cases that demonstrate the utility and The elemental and isotopic data for these studies were from ostracods identified as 428 Candona sp., by R. Forester (see discussion in Curtis and Hodell, 1993), which is an ostracod 429 species closely related to Candona rawsoni. See table 2 for detailed model input parameters, 430 including the partition coefficients for Candona rawsoni as determined by Engstrom and Nelson 431 (1991). 432 The modeled temperature and lake volume changes for the Lake Miragoane record are 433 largely in accordance with the expected response (see interpretation in Curtis and Hodell, 1993, 434 their fig. 8). These authors suggested based on the isotopic and elemental data that lake levels 435 Ibarra and Chamberlain -21 were highest and climate was "mesic" (wetter and hotter) during the early Holocene, peaking 436 between 7,000 and 4,000 years BP (Hodell and others,1991; Curtis and Hodell, 1993). In the late 437 Holocene they suggest that between ~4,000 and ~2,500 years BP there was a two-step increase in 438 aridity towards lower lake level and increased salinity, based primarily on the oxygen isotopes 439 (Hodell and others, 1991), followed by wetter conditions beginning at ~1,000 year BP based on 440 additional evidence from the trace element measurements (Curtis and Hodell, 1993). 441 Our model agrees with this interpretation and shows that modeled temperature and lake 442 volume covary. Lake volume and temperature increase to a peak at ~6300 years BP, decline to a 443 minimum at ~2000 years BP, and increase again nearing present day. The deglacial warming 444 (~10,000 to 7,000 years BP), and increasing aridity (from ~4,000 to 1,000 years BP) are captured 445 by the predicted temperature fluctuations. Driven by the covariation in the input datasets, the 446 modeled temperature is substantially different than that expected from just using the 0%. However, we point out that due to the large spread in the measured Sr/Ca and Mg/Ca, the 496 modeled volume changes are not particularly robust. 497

LIMITATIONS AND GUIDELINES 499
The modeling framework presented in this paper represents a new approach to 500 quantitatively constraining the hydrologic balance and temperature fluctuations as recorded by 501 the trace elements and stable isotopes of biogenic lacustrine carbonate. Inherently the modeling of paleoclimate records requires simplifying assumptions. As such, we provide several 503 guidelines for assessing the potential application of HyBIM to a lake system: 504 (1) Lake systems with highly variable (inorganic) chemical sedimentation, typically driven 505 by order of magnitude changes in lake level, cannot be modeled using HyBIM. A 506 fundamental assumption of the evaporation partial derivatives is that we assume mass 507 conservation of the trace elements within each time step. In addition, we assume that 508 within each time step lake volume does not change by more than 10% (see derivation of 509 partial derivatives). Thus, HyBIM solutions should be limited to "balanced-filled" (Caroll 510 and Bohacs, 1999) lake systems with biogenic carbonate from periods of similar 511 lithologic deposition. 512 (2) Many large lakes break up into smaller lakes during desiccation due to complex 513 hypsometric basin geometries and watersheds (for example Quaternary Lake Lahontan; 514 Reheis and others, 2014). While this is not readily assessed from sediment cores alone, HyBIM to a high-resolution historical record from a small, constrained terminal lake system will 554 provide the best independent validation of this modeling methodology. Finally, the Monte Carlo approach presented in this paper only includes uncertainty in the 587 input parameters and input dataset. However, age-models in sedimentary settings, and many 588 terrestrial paleoclimate archives, including lakes, can be highly uncertain due to sedimentation 589 rate changes, hiatuses and a lack of dateable material (for example, Huybers and Wunsch, 2004;