USING δ2H AND δ18O TO DETERMINE THE FLOWPATHS AND TIMESCALES OF WATER AT THE SUSQUEHANNA SHALE HILLS CRITICAL ZONE OBSERVATORY
ABSTRACT
A stable isotope sampling network was implemented at the Susquehanna Shale Hills Critical Zone Observatory. The objective was to determine the δ2H and δ18O signature in the catchment pools to determine the flowpaths and timescales of the hydrologic system. The stable isotope network covers all phases of the hydrologic cycle, including precipitation sampled adaptively during precipitation events with an Eigenbrodt NSA-181/S wet-only collector (sixhour samples), soil water sampled weekly along four transects with 80 suction-cup lysimeters, groundwater sampled daily at two wells with ISCO automatic samplers and bi-weekly at 16 wells and stream water sampled daily with an ISCO automatic sampler. The comprehensive sampling of the network was possible because of the DLT-100 liquid water stable isotope analyzer from Los Gatos Research, with a reproducibility of ± 0.2%0 for δ18O, ± 1.0%0 for δD and the capability to run 30 samples per day. The δ2H and δ18O data showed the dominance of cold season infiltration and recharge, with recharge specifically occurring over the period of late September – May. The δ2H and δ18O record also showed that groundwater regularly flushed the deep soil water, and that groundwater is the major component of streamflow. Preferential flowpaths in the soil during the cool or non-growing season was identified and is related to stream stormflow. A piecewise constant model for flow, tracer concentration and age was based on the work of Duffy and Cusumano (1998) and Duffy (2010), and was unique in that it solved for transient flow. The finding of the age model was that the mean age of the water in the catchment ranged between 4.5 – 9 months. The oldest ages occurred during the summer drought and the youngest ages occurred during times of maximum recharge over the winter. This research was performed as part of the NSF-funded Critical Zone Observatory and the importance of this effort multi-investigator effort was essential to the success of this research.
TABLE OF CONTENTS
LIST OF FIGURES……………………………………………………………………………………………………v
LIST OF TABLES…………………………………………………………………………………………………….vi
ACKNOWLEDGEMENTS………………………………………………………………………………………..vii
Chapter 1 Introduction ………………………………………………………………………………………………1
Stable Isotope Overview ……………………………………………………………………………………..2 δ2H and δ18O in Catchment Hydrology………………………………………………………………….4
Age Modeling Theory …………………………………………………………………………………………9
Chapter 2 Site Description and Sampling Methods……………………………………………………….13
Sampling Methods………………………………………………………………………………………………13
SSHCZO Conceptual Flow Model………………………………………………………………………..16
Chapter 3 δ2H and δ18O Results ………………………………………………………………………………….19
Meteoric Water Line …………………………………………………………………………………………..19
Precipitation ………………………………………………………………………………………………………22
Soil Water………………………………………………………………………………………………………….26
Groundwater………………………………………………………………………………………………………29
Stream Water……………………………………………………………………………………………………..32
Chapter 4 δ2H and δ18O Discussion ……………………………………………………………………………..35
Chapter 5 Age Model………………………………………………………………………………………………..39
Derivation………………………………………………………………………………………………………….39
Age Model Results and Discussion……………………………………………………………………….44
Chapter 6 Conclusion………………………………………………………………………………………………..52
References………………………………………………………………………………………………………………..54
Appendix δ2H and δ18O Data and Age Model Output …………………………………………………….59
Chapter 1
Introduction
The Susquehanna Shale Hills Critical Zone Observatory (SSHCZO) was designed for the ultimate goal of determining the formation, evolution and function of regolith. The functioning of regolith affects the hydrologic flowpaths and timescales of the water in a catchment (Anderson et al., 2008; Brantley et al., 2007). Hydrologic studies have occurred at Shale Hills since the 1970s focusing on streamflow generation (Lynch, 1976; Duffy, 1996) soil moisture and preferential flow (Lin, 2006; Lin et al., 2006; Lin and Zhou, 2008; Graham and Lin, 2011) and solute transport (Jin et al., 2010; Jin et al., 2011; Kuntz, 2010). These studies have found that the streamflow largely depends on the antecendent soil moisture and groundwater, which is why the stream does not flow over the summer growing season due to low soil moisture and a low water table (Lynch, 1976). Once the soil moisture deficit is satisfied than groundwater recharge can occur over the cold season or non-growing season (Deines et al., 1990; O’Driscoll et al., 2005).
Subsurface water is transported through several preferential flowpaths at Shale Hills (Lin, 2006;
Lin et al., 2006; Lin and Zhou, 2008), and the groundwater regularly flushes the deep soil water (Lynch and Corbett, 1989; Duffy, 1996). A conceptual hydrologic model can be pieced together from these studies but no study has looked at the hydrology as a whole. The goal of this research is to identify the composition of the hydrologic pools and the flowpaths and timescales through the catchment and to evaluate how these pools interact and mix to form runoff within the watershed. To accomplish this an isotopic sampling network was developed for the collection of δ2H and δ18O samples. δ2H and δ18O were chosen because they are conservative tracers at low temperatures (Hoefs, 2009) and they compose the water molecule. The age of the hydrologic system was determined by modeling the flow and isotopic concentration of the stream. The model used was created by Duffy (2010), and was chosen because it is the first of this type of model to solve for transient flow.
Stable Isotope Overview
Isotopes are measured as a ratio of the heavy (rare) to light (common) isotope, e.g. R = 2H/1H or 18O/16O. Measuring the exact ratio of a sample is inaccurate because the heavy isotopes are not abundant, instead it is orders of magnitude more accurate to relate the ratio of a sample to a known ratio of standard. The standard used for water is Vienna Standard Mean Ocean Water
(VSMOW), and the absolute abundance of deuterium (2H or D) in VSMOW is 2H/1H = 155.95 x 10-6 (Dewit et al., 1980), and 18O is 18O/16O = (2005.2 ±0.45)x 10-6 (Baertschi, 1976). Isotope values are reported in delta (δ) notation
Rsample −Rstandard
δsample =[ ] ×1000
Rstandard
and the values are on a permil (‰) scale. A positive δ is enriched in the heavy isotope compared to the standard and a negative δ is depleted in the heavy isotope compared to the standard.
The abundance of stable isotopes varies in nature due to fractionation. Fractionation occurs because of differences in mass between isotopes, which causes the isotopes to have different physiochemical properties. The physical properties of an isotope affect its vibrational motions only, therefore causing differences in zero-point energies of isotopes (Hoefs, 2009). Therefore it is harder to break a bond between two 2H isotopes then it is between two 1H isotopes, because the bond of the light isotope is weaker. Fractionation is significant at low temperatures and disappears at higher temperatures because of the amount of energy available. Fractionation is measured using the fractionation factor, α, or epsilon, ε.
RA 1000+δA αA−B = ≈
RB 1000+δB
εA−B =(αA−B −1)×1000
Fractionation is important for light isotopes because of larger mass differences. For example, using the atomic weights of Sharp (2007), 2H is 99.8% heavier then 1H, 18O is 12.5% heavier then 16O, but a heavier isotope like 238U is only 1.27% heavier then 235U.
Fractionation occurs as either an equilibrium or kinetic reaction (Hoefs, 2009). Equilibrium fractionation occurs when the compound is in a closed environment and the products and reactant are able to reach isotopic equilibrium. The backwards and forwards reaction rates are equal, and the δ value in each compound reaches a constant value. In equilibrium fractionation it is generally understood that the heavy isotopes will preferentially accumulate in the denser state (eg. solid>liquid>vapor).
Kinetic fractionation is an irreversible process, meaning that the backwards and forwards reaction rates are not equal. This can be due to the products being isolated from the reactants or the products being carried away from the reactants, like in the case of evaporation. The amount of fractionation depends on the ratio of masses of the isotopes in the products and reactants and their bond strengths. It is generally understood that during kinetic fractionation the light isotopes accumulate in the products because it takes less energy to break their bonds.
The most important fractionations with respect to water are related to phase changes, due to differences in vapor pressure and freezing point (Friedman et al., 1964). Water molecules containing 18O and 2H cause the water to have a lower vapor pressure and freezing point, therefore preferentially accumulating in the more dense phases.
When studying the stable isotopes of water it is only necessary to consider water
molecules of the form 1H216O, 1H2H16O, and 1H218O. This is because these forms occur in concentrations that are orders of magnitude larger then the concentrations of the other six forms of the water molecule (Friedman et al., 1964).
δ2H and δ18O in Catchment Hydrology
δ2H and δ18O isotopes are ideal in catchment studies because the isotopes fractionate and mix as they move between the different pools in the catchment. This has led to δ2H and δ18O being used for hydrograph separations (Skalsh et al., 1976; Sklash and Farvolden, 1979; Rice and Hornberger, 1998), estimation of soil water movement and evaporation (Allison and Barnes, 1983; Barnes and Allison, 1988; Gazis and Feng; 2004), estimation of groundwater composition and recharge (Gat 1971; Davisson and Criss, 1993; Criss and Davisson, 1996; Darling and Bath, 1988; Winograd et al., 1998) and the determination of preferential flowpaths and the old water paradox (McDonnel, 1990; Leaney et al., 1993; Kirchner, 2003; Vogel et al., 2010). It is necessary to understand how δ2H and δ18O function in catchments so that these isotopes can be used for these purposes.
Precipitation is the only input to the SSHCZO, so the δ2H and δ18O signature of precipitation serves as the starting point for the δ2H and δ18O signature that will be acquired by the soil water, groundwater and stream water. The main factors controlling the δ2H and δ18O of a precipitation event is the source of water, temperature during condensation, fraction of original water remaining, amount of recycled (evapotranspired) water, pathway of the event and isotopic exchange of water droplets and water vapor.
The source water region is important because of the temperature and relative humidity at the time of evaporation. A higher temperature will decrease the fractionation factor while the relative humidity determines the amount of evaporation and exchange between the water vapor and liquid water (Merlivat and Jouzel, 1979). Evaporation is a kinetic fractionation process, and the water vapor is always depleted in δ2H and δ18O relative to the liquid water (Dansgaard, 1964). The water vapor receives more δ2H than δ18O though because it weighs less and therefore takes slightly less energy to break the bond.
Condensation is an equilibrium fractionation process where the amount of fractionation is controlled by the air temperature (Dansgaard, 1964). Because of this the δ2H and δ18O of precipitation correlates with temperature on a global scale (Dansgaard, 1964; Siegenthaler and Oeschger, 1980; Lawrence, 1980; Rozanski et al., 1992). Since condensation is an equilibrium process, theoretically the initial condensate will have the same composition as the source water, but this is not the case due to kinetic evaporation (Dansgaard, 1964).
As condensation continues on a finite volume of water vapor, the heavy isotopes become depleted in the water vapor because they are more stable in the liquid water. This leads to four ‘effects’, which were discovered by Dansgaard (1964) and Friedman et al. (1964). They are the amount effect, continental effect, latitude effect and altitude effect. The premise behind all of the effects is that lower temperatures and a finite volume of water vapor lead to increased fractionation, which occurs as a spatial or temporal gradient as rainout occurs, or an air mass moves over the continent, etc.
Precipitation isotopic values are also influenced by the amount of water evapotranspired back to the air mass. Concentration gradients calculated for the Amazon basin show a reduced ‘continental effect’ when compared to other regions (Araguas-Araguas et al., 2000). Unlike evaporation transpiration does not fractionate water so the water vapor returning to the air mass is comparatively enriched.
During decent liquid precipitation is subject to isotopic alteration by evaporation and isotopic exchange with atmospheric moisture. Evaporation occurs at the beginning of a precipitation event, because the atmosphere underneath the air mass is normally not saturated with water, and after saturation of the underlying air isotopic exchange occurs. The amount of time or distance it takes for a water drop to isotopically equilibrate with the atmosphere is positively related to the drop radius and the atmospheric temperature (Friedman et al., 1962). It should be noted that during heavy rain, it is the liquid component that dominates the exchange process causing the atmospheric water vapor to have the same δ2H and δ18O as the liquid (Dansgaard, 1964). Isotopic exchange also affects other open water bodies. The rate of exchange is controlled by the volume and surface area of the water bodies (Ingraham and Criss, 1993). In a study using beakers of different volume and surface area Ingraham and Criss found differences of up to 9% of the original isotopic composition after five days.
The pathway of a precipitation event determines the conditions under which condensation will occur. In a study of precipitation δ2H and δ18O and the pathway of the events at Mohonk, NY, Lawrence et. al (1982) found a distinct difference among the isotopic signatures and the corresponding path of the precipitation event. Their finding was that as the center of the storm is displaced seaward the δ2H of the precipitation decreases. The reason for the decrease is because the frontal surface creating these events is at a higher altitude when the storms are seaward.
Isotopic values of precipitation also vary during a storm. Pionke and DeWalle (1992) collected precipitation samples over 21 minute intervals for 33 storms in central Pennsylvania and found that short storms showed little change in the δ2H and δ18O, but long storms showed great variability that seemed to be related to precipitation intensity. Differences in δ18O of 15‰ were found for an 11-hour storm that dropped 31 mm of precipitation. This intra-storm variability is important because a specific part of the storm may dominate the hydrologic response of the catchment.
Precipitation δ2H and δ18O is variable because it is a global scale process, but this variation is attenuated in the other hydrologic pools. Generally the soil water, groundwater and stream water pools are isotopically affected by fractionation from evaporation and freeze/thaw and mixing with water bodies that have a different δ2H and δ18O signature. Evaporation enriches the δ2H and δ18O signature of water, freezing depletes the δ2H and δ18O signature in liquid water and melt water gradually enriches the δ2H and δ18O signature. Freezing and melting are both kinetic fractionation processes but they leave no permanent change on the δ2H and δ18O composition as long as all of the water is frozen or melted.
There are many ways that δ2H and δ18O signatures can be used to understand the hydrologic functioning of a catchment. Hydrograph separations are possible because of the difference between precipitation δ2H and δ18O and that of the soil water, groundwater and stream water. Sklash et al. (1976) and Sklash and Farvolden (1979) used a simple linear two-component mixing model to determine the relative contributions to streamflow from event and pre-event water. Pre-event water was any water found in the catchment before the precipitation event while the event water was the precipitation. These studies were the first to show that pre-event water, the groundwater and soil water, made up the bulk of the streamflow during a precipitation event. More recent studies have used three-component and even five-component models, by using additional tracers, to determine the relative contributions of soil water, groundwater and precipitation (Rice and Hornberger, 1998; Brown et al., 1999; Uhlenbrook and Hoeg, 2003;
Sayama and McDonnell 2009).
The finding that pre-event water contributes the majority of water to the stream during a storm has led to the old water paradox (Kirchner, 2003). The paradox is that during a precipitation event “old” subsurface water is quickly mobilized and forms the bulk of the streamflow, rather then the precipitation falling on the catchment. So far scientists have tried to explain why this water is mobilized and not the event water, and McDonnell (1990) even used δ2H for his study. McDonnell’s hypothesis was that the old water was quickly transported to the stream through macropores. The reason that the bulk of the water was old water was because the soil was almost completely saturated. Therefore he argued that it only took a small amount of event water to saturate the soil and initiate macropore flow. The old water paradox still puzzles scientists, despite many attempts to explain it.
δ2H and δ18O have proven useful in estimating subsurface evaporation and flowpaths. Allison and Barnes (1983) modeled soil water evaporation using the enrichment of δ2H and showed that evaporation typically occurs in the top one meter of soils. It is possible to trace the movement of soil water due to isotopically distinct precipitation events. The infiltrated precipitation serves as an isotopic tag that can be followed and used to determine infiltration and recharge rates and soil water flowpaths (Barnes and Allison, 1988; Bengtsson et al., 1987; Darling and Bath, 1988; Gazis and Feng, 2004). δ2H and δ18O have also been used to understand the functioning of preferential flowpaths and mobile and immobile water (McDonnell, 1991; Leaney et al. 1993; Vogel et al., 2010; Brooks et al., 2010). In particular, Leaney et al. (1993) found that water travelling through preferential flowpaths was bypassing parts of the soil column and contributing a large portion of water to streamflow. Brooks et al. (2010) have also found that mobile and immobile water have different δ2H and δ18O signatures, which could change the way the unsaturated zone is conceptualized.
δ2H and δ18O have also been used to investigate groundwater flowpaths and the timing and constituents of recharge water (Gat, 1971). Davisson and Criss (1993) were able to show that the composition of recharge changed over the course of the year depending on the relative height of the groundwater. Winograd et al. (1998) were also able to determine the make up of groundwater recharge by determining the precipitation δ2H and δ18O seasonal signatures. δ2H and δ18O in hydrologic studies are normally plotted on a meteoric water line plot (mwl). Craig (1961) was the first to plot δ2H vs δ18O for water samples collected around the world and find that δ2H and δ18O are linearly related. The linear relationship formed by global precipitation samples is called the global meteoric water line and has the relationship
δ2H=8δ18O+10.
Local precipitation samples will have a slightly different linear relationship, and they are called local meteoric water lines. The significance of the lmwl is that fractionated water samples do not plot along the lmwl. Evaporated samples plot below the lmwl because slightly more δ2H is evaporated relative to δ18O (Hoefs, 2009). Snowmelt plots above the lmwl because slightly more δ2H is melted relative to δ18O (Clark and Fritz, 1997). Therefore the mwl plot is essential to the interpretation of δ2H and δ18O data from catchment waters.
Age Modeling Theory
The theory of age and residence time modeling for watersheds was developed through the early work of chemical engineers modeling tank reactors (Danckwerts, 1953; Nauman, 1969, 1973). Eriksson (1971) and Bolin and Rodhe (1973) gave concise summaries of the reservoir theory, which brought the theory of age modeling clearly into the field of hydrology. The premise behind the reservoir theory is that an element entering a well-mixed reservoir is characterized by τ, the time that has elapsed since it entered the reservoir. The elements are then arranged in a cumulative function M(τ), or F(τ) for the flux leaving the reservoir, which is the total mass of elements in the reservoir that has spent a time equal to or less then τ in the reservoir. M(τ) or F(τ) can be thought of as the residence time distribution (rtd), which gives the distribution of residence times for all of the elements in the reservoir Then the steady-state mean residence time or turnover time is defined as τ0 = M0/F0, the total mass divided by the total flux. The average age of elements in the reservoir is defined as
1 ∞
τa = M0 ∫0 τdM(τ).
The rtd of a reservoir is important because a mean residence time can be achieved many ways, so the rtd shows how the reservoir is unique. The reservoir theory only applies to the steady state flow case, as do most age models.
Maloszewskie and Zuber (1982) introduced several lumped parameter models for the determination of the residence time from isotope data. The convolution integral was used along with a weighting function that specified how the isotope moved through the system:
t
Cout Cin (τ)exp[−λ(t −τ)]g(t −τ)dτ
where Cout and Cin are the output and input concentrations, τ is the transit time of the isotope, λ is the half-life for a radioactive isotope, and g is the weighting function. There are several weighting functions that can be used including piston flow, exponential flow, combined exponential and piston flow, and dispersive flow. The correct weighting function is chosen by matching the modeled output concentration to the measured output concentration. The transit time or age of the tracer for this model was defined as
dt
. dt
The convolution integral has been widely used even though there are several drawbacks. First is that it is difficult to accurately define the input concentration, and the input is supposed to be instantaneous. Bergmann et al. (1986) developed an equation for the input function that uses α, an infiltration coefficient, to better define the input concentration. A second problem is that this model assumes steady state flow, which is not found in nature. Werner and Kadlec (1996) show that the output concentration from a system depends on when the tracer entered the system relative to the flow through the system. To get around this problem they use a flow-weighted time, which weights the instantaneous volume by the total volume that leaves the system during the period of interest. Zuber (1986) also got around the steady state assumption by rewriting the lumped parameter model equations using the tracer mass flux. He notes though that this can only work for catchments with short residence times because a longer input and output concentration time series is required. The general lumped parameter model can also not deal with mobile and immobile flow and preferential flow. The advantage of the lumped paramter model of Maloszewski and Zuber (1982) is that it gives the rtd of the water rather then just a mean residence time.
Goode (1996) was able to directly simulate the age of groundwater by developing an adevective – dispersive equation that conserves the conceptual “age mass”. Age mass is defined as the mean age of a parcel of water multiplied by the mass of the water parcel. Therefore when two water bodies mix the age is a mass-weighted average of the two water bodies. More recently Delhez et al. (1999) and Gourgue et al. (2006) have shown that by using age mass it is possible to analyze the moments of the concentration equation without actually needing to define the exact form of the equation. This allowed Duffy (2010) to develop a coupled dynamical system of differential equations that include transient flow and a theoretical approach to tracer age that does not require the residence time functional form, only the first 2 moments. The system of equations was formed from general equations of flow and concentration for a volume averaged system. The total fluid volume of the system satisfied a balance equation:
dV dt
= Qin −Qout
and a tracer mass balance equation:
dVC dt
= QinCin −QoutCout +VΓc
where Γc represents the internal sources and sinks for the tracer. The final dynamical set of equations that came from the balance equations are:
dV dt
= Qin −Qout
dC = Qin (Cin −Cout )+Γcdt V
dα= Qin
Cout − α+Γα dt V
where α is the age concentration. The tracer mean age is defined as
α(t)
A(t)= .
C(t)
Duffy (2010) extended the dynamical system to include mobile and immobile flow, which have been shown to be isotopically distinct (Brooks et al., 2010).