In cvoter/isoH2Obudget: Lake Water Budget Using Stable Isotopes
# Libraries
library(devtools)
library(CSLSiso)
library(CSLSdata)
library(reshape2) # for melting
library(ggplot2) # for plotting
library(cowplot) # for subplots
library(extrafont) # for fonts
library(dplyr) # for %>% pipes
library(lubridate)
library(knitr)
library(stringr)
library(scales) # for ::
# Parameters
text_size <- 12
lakes <- c("Pleasant", "Long", "Plainfield")
# Speed up processing by getting combinded inputs and outputs once here.
inputs <- NULL
outputs <- NULL
for (lake in lakes){
# Inputs
these_inputs <- runall_csls_inputs(lake)
these_inputs$lake <- lake
inputs <- rbind(inputs, these_inputs)
# Outputs
these_outputs <- runall_csls_budget(lake)
these_outputs$lake <- lake
outputs <- rbind(outputs, these_outputs)
}
inputs$lake <- factor(inputs$lake, levels = lakes)
outputs$lake <- factor(outputs$lake, levels = lakes)
Overview {#top}
This water budget model uses 1) observed changes in lake levels, 2) observed
precipitation, 3) calculated lake evaporation, and 4) stable isotope
measurements of precipitation, evaporation, lake water, and inflowing
groundwater to calculate monthly groundwater flow into and out of a lake.
The lakes in the Central Sands Lakes Study (Pleasant Lake, Long Lake, and
Plainfield Lake) are all seepage lakes, a type of lake that does not have
inflowing streams, outflowing streams, or appreciable surface runoff. The
water budget for a seepage lake is as follows:
\begin{align}
\tag{1}
\frac{dV}{dt} = P - E + GW_{in} - GW_{out}
\end{align}
where $\frac{dV}{dt}$ is the change in lake volume (m^3^), $P$ is the
precipitation (m^3^), $E$ is the lake evaporation (m^3^), $GW_{in}$ is the
groundwater flowing in to the lake (m^3^), and $GW_{out}$ is the groundwater
flowing out of the lake.
Change in Lake Volume is directly measured using a water level
logger at each lake. Daily lake levels elevations are converted into daily lake
volumes using lake elevation-volume relationships developed from lake bathymetry
maps.
Precipitation is directly measured at a nearby weather station.
It is converted to a volume by multiplying the total precipitation (mm) by the
average lake area during the month of interest. Average lake area is determined
using daily lake level elevation from water level loggers and the lake
elevation-area relationships developed from bathymetry maps.
Lake Evaporation is calculated based on the lake energy-budget
approach outlined in McJannet et al. (2008) as described in McMahon et al.
(2013). This approach uses:
- Measured incoming solar radiation, wind speed, relative humidity, and air
temperature from a nearby weather station
- Lake surface temperature from continuous temperature loggers
- Area and depth of the lake from daily lake stage measurements and
elevation-area relationships
Evaporation as a depth (Eq. 2) is summarized at a monthly timestep and converted
to a volume using the average lake area during that month, as in the
precipitation calculations.
\begin{align}
\tag{2}
E = \frac{\Delta_{w}(R_{n} - G_{w}) +
\frac{606024\rho_{a}c_{a}(es_{w} - es_{a})}{r_{a}}}{\lambda*(\Delta_{w} + \gamma)}
\end{align}
where $E$ is lake evaporation (mm/day), $\Delta_{w}$ is the slope of the
saturation water vapour curve at water temperature (kPa/deg C), $R_{n}$ is the
net radiation (MJ/m^2^/day), $G_{w}$ is the change in heat storage in the water
body (MJ/m^2/day), $\rho_{a}$ is the density of air (kg/m^3^), $c_{a}$ is the
specific heat of air (MJ/kg/K), $es_{w}$ is the saturation vapour pressure at
water temperature (kPa), $es_{a}$ is the saturation vapour pressure of the air
(kPa), $r_{a}$ is the aerodynamic resistence (s/m), $\lambda$ is the latent heat
of vapourisation (MJ/kg), and $\gamma$ is the psychormetric constant (kPa/deg
C).
Groundwater inflow is calculated using monthly changes in lake
volume, precipitation, lake evaporation, and stable isotope measurements based
on the approach in Krabbenhoft et al. (1990) and modified based on Mook (2000) to
account for non-steady state conditions in lake levels:
\begin{align}
\tag{3}
GW_{in} = \frac{P(\delta^{18}O_{L} - \delta^{18}O_{P}) + E(\delta^{18}O_{E} -
\delta^{18}O_{L}) + V*d\delta^{18}O_{L}/dt}{\delta^{18}O_{GW_{in}} - \delta^{18}O_{L}}
\end{align}
where $GW_{in}$ is the groundwater inflow (m^3^), $P$ is precipitation (m^3^), $E$
is evaporation (m^3^), V is the lake volume (m^3^) and $\delta^{18}O_{x}$ is the
relative isotopic composition in units of per mil relative to a known standard.
$\delta^{18}O_{L}$ is measured at the lake, $\delta^{18}O_{P}$ is measured from
precipitation, $\delta^{18}O_{E}$ is estimated for evaporation, and
$\delta^{18}O_{GW_{in}}$ is measured at inflowing groundwater wells.
Groundwater outflow can then be calculated as the only remaining
term in the water balance:
\begin{align}
\tag{4}
GW_{out} = GW_{in} + P - E - \frac{dV}{dt}
\end{align}
Detailed Approach
Input data for the Central Sands Lake Study sites (Pleasant Lake, Long Lake, and
Plainfield Lake) are saved in the CSLSdata R data package. The CSLSdata
package includes a data-raw directory where raw data is stored together with
the scripts used to clean raw data and save them as Rda files. Detailed
information on the source of each dataset is captured in the CSLSdata manual.
This is accessible as a pdf or via RStudio (help(package = CSLSdata)).
Below, we describe all additional assumptions and calculations used to transform
CSLSdata Rda files from their form in the CSLSdata package into the monthly
inputs required for monthly water balance calculations.
Change in Lake Volume {#dV}
We obtain daily lake elevation data from water level loggers at USGS gages
at all three lakes in the Central Sands Lakes Study.
plot_lake_level(text_size = text_size)
We also have the elevation-volume relationships for each lake from lake
bathymetry data (calculated in ArcGIS).
plot_elev_area_vol(ydata = "vol_m3",
ylabel = expression(Lake~Volume~(m^{3})),
text_size = text_size)
We convert these daily lake level elevations to daily lake volumes using these elevation-volume relationships.
plot_lake_level(ydata = "vol_m3",
ylabel = expression(Lake~Volume~(m^{3})),
yscales = scales::scientific,
text_size = text_size)
The change in lake volume for a given month is then the difference
between the lake volume on the last day of the month and the lake volume on the
first day of the month.
plot_input(df = inputs,
ydata = "dV_m3",
ylabel = expression(Change~Lake~Volume~(m^{3})),
geom_type = "col",
text_size = text_size)
Precipitation {#precip}
We obtain hourly precipitation depth (mm) from the Hancock Agricultural
Research Station (station
id: hck). The station is located in Hancock, WI (location: 44.1188, -89.533,
elevation: 241m), approximately 8 miles from Plainfield Lake, 8.5 miles from
Long Lake, and 14.5 miles from Pleasant Lake.
weather <- CSLSdata::weather
monthly_weather <- weather %>%
group_by(date = floor_date(.data$date, unit = "month")) %>%
summarise(P = sum(.data$P))
ggplot(monthly_weather, aes(x = date, y = P)) +
geom_col() +
scale_x_datetime(date_breaks = "3 months",
date_labels = "%b '%y") +
scale_y_continuous(labels = scales::scientific) +
labs(x = "",
y = "Monthly Precipitation (mm)") +
theme_bw() +
theme(text = element_text(family = "Segoe UI Semilight",
size = text_size),
plot.title = element_text(hjust = 0.5),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
To convert precipitation as a depth to precipitation as a volume, We use the
daily lake elevations at the USGS gages at each lake.
plot_lake_level(text_size = text_size)
We also use the elevation-area relationships for each lake from lake
bathymetry data (calculated in ArcGIS).
plot_elev_area_vol(ydata = "area_m2",
ylabel = expression(Lake~Area~(m^{2})),
text_size = text_size)
We convert daily lake level elevations to daily lake areas using these
elevation-area relationships.
plot_lake_level(ydata = "area_m2",
ylabel = expression(Lake~Area~(m^{2})),
yscales = scales::scientific,
text_size = text_size)
Then, we calculate the mean monthly lake area.
plot_input(df = inputs,
ydata = "mean_area_m2",
ylabel = expression(Lake~Area~(m^{2})),
geom_type = "point",
facet_scales = "free_y",
text_size = text_size)
By multiplying monthly precipitation as a depth (mm) by the mean monthly lake
area (m^2^), we obtain monthly precipitation as a volume (m^3^) for each lake.
plot_input(df = inputs,
ydata = "P_m3",
ylabel = expression(Precipitation~(m^{3})),
geom_type = "col",
text_size = text_size)
Lake Evaporation {#evap}
All lake evaporation calculations are contained in the CSLSevap package; more
details on this energy-balance approach can be obtained by searching the
documentation for that package. Broadly, the energy-budget lake evaporation
equations require information on:
- Lake location (latitude, longitude, elevation, latitude of time zone)
- Lake size (area, depth)
- Lake surface temperature (daily values, taken from surface-most continuous
temperature logger)
- Meteorological parameters (daily relative humidity, air temperature, incoming solar
radiation, wind speed from nearby weather station)
- Heat absorption and reflection parameters (e.g., albedo of water, density and
specific heat of air)
These lake evaporation equations yeild daily lake evaporation (mm/day), which is
converted to monthly lake evaporation volume by summing over each month and
multiplying by the average monthly lake area, as in the precipitation
calculations.
plot_input(df = inputs,
ydata = "E_m3",
ylabel = expression(Evaporation~(m^{3})),
geom_type = "col",
text_size = text_size)
Groundwater Inflow {#GWin}
As indicated in Eq. 3, calculation of groundwater inflow requires:
-
Monthly precipitation, evaporation, and mean lake volume. These
are calculated per the previous sections.
-
Monthly groundwater stable isotope measurements are derived from regular
CSLS sampling at four groundwater wells near each lake. These measurements were
taken approximately every other month. Upgradient wells representing inflowing
groundwater are defined as wells where the mean daily difference between
groundwater levels and lake levels in the 30-day period prior to sampling is
greater than 1 cm (i.e., the precision of the water level logger). All
upgradient well measurements from the same sampling event for a lake are
averaged, then linearly interpolated to estimate values on the last day of each
month.
-
Monthly lake stable isotope measurements are derived from regular CSLS
sampling at each lake. These measurements were taken less frequently than
groundwater measurements; approximately every other month during the ice-free
period. All lake measurements are linearly interpolated to estimate values on
the last day of each month during the ice-free period.
-
Monthly precipitation stable isotope measurements are derived from
regular CSLS sampling at two precipitation collectors located at Hancock
Agricultural Research Station. These measurements were taken approximately every
month during the snow-free season. They are supplemented by monthly measurements
taken by Maribeth Kniffin at the same location from May 2016 to April 2017.
Since 1) we do not have access to the exact days that Maribeth Kniffin collected
her measurements and 2) precipitation samples are an aggregate sample of all
precipitation collected between sample dates (i.e., they are not a snapshot in
time), we do not perform any linear interpolation on these samples. Instead, we
assume all precipitation measurements represent precipitation on the last day of
the month in which they were collected. Maribeth Kniffin's data is only used for
months when we do not have CSLS-collected data.
lakes <- c("Pleasant", "Long", "Plainfield")
monthly_isotopes <- NULL
for (lake in lakes){
isotopes <- CSLSdata::isotopes[[lake]]
dictionary <- CSLSdata::dictionary[[lake]]
lake_levels <- CSLSdata::lake_levels[[lake]]
gw_levels <- CSLSdata::gw_levels[[lake]]
# Measured values
tmp1 <- iso_site_type(isotopes, dictionary, lake_levels, gw_levels)
tmp1 <- tmp1 %>%
group_by(floor_date(.data$date, unit = "month"),
variable = .data$site_type) %>%
filter(is.na(.data$site_type) == FALSE,
.data$site_type != "",
.data$site_type != "downgradient",
.data$site_type != "deep") %>%
summarise(date = floor_date(mean(.data$date), unit = "day"),
value = mean(.data$d18O)) %>%
ungroup() %>%
select(.data$date, .data$variable, .data$value)
tmp1$variable <- tmp1$variable %>%
str_replace("upgradient","Groundwater") %>%
str_replace("lake","Lake") %>%
str_replace("precipitation","Precipitation")
tmp1$lake <- lake
tmp1$measurement <- "actual"
# monthly values
tmp2 <- summarise_isotopes(isotopes, dictionary,
find_timeseries(isotopes),
lake_levels, gw_levels) %>%
select(date, d18O_GWin, d18O_lake, d18O_pcpn) %>%
melt(id.vars = "date")
tmp2$variable <- tmp2$variable %>%
str_replace("d18O_GWin","Groundwater") %>%
str_replace("d18O_lake","Lake") %>%
str_replace("d18O_pcpn","Precipitation")
tmp2$lake <- lake
tmp2$measurement <- "monthly"
tmp3 <- tmp1 %>% filter(variable == "Precipitation")
tmp3 <- fill_timeseries_gaps(tmp3, find_timeseries(isotopes))
kniffin_isotopes <- CSLSdata::kniffin$isotopes
for (i in 1:nrow(tmp3)) {
this_month <- month(tmp3$date[i])
this_i <- which(month(tmp3$date) == this_month)
tmp3$value[i] <- mean(tmp3$value[this_i], na.rm = TRUE)
tmp3$measurement[i] <- "actual"
tmp3$value[is.nan(tmp3$value)] <- NA
this_kniffin <- kniffin_isotopes %>%
filter(month(.data$date) == this_month,
.data$site_id == "PRECIP")
if (nrow(this_kniffin) == 0) {
kniffin_d18O <- NA
} else {
kniffin_d18O <- mean(this_kniffin$d18O, na.rm = TRUE)
}
if (is.na(tmp3$value[i])){
tmp3$value[i] <- kniffin_d18O
tmp3$measurement[i] <- "kniffin"
}
}
tmp3$lake <- lake
tmp3$variable <- "Precipitation"
monthly_isotopes <- rbind(monthly_isotopes, tmp1, tmp2, tmp3)
}
monthly_isotopes$lake <- factor(monthly_isotopes$lake, levels = lakes)
ggplot(monthly_isotopes %>% arrange(measurement),
aes(x = date, y = value, color = measurement, shape = measurement)) +
geom_point() +
scale_x_datetime(date_breaks = "4 months",
date_labels = "%b '%y") +
facet_grid(variable~lake) +
labs(x = "", y = expression(delta^18*'O'), title = "") +
scale_color_manual(name = "",
values = c("red", "purple", "grey70"),
labels = c("Measured by CSLS",
"Measured by M. Kniffin",
"Monthly Interpolation")) +
scale_shape_manual(name = "",
values = c(19, 15, 17),
labels = c("Measured by CSLS",
"Measured by M. Kniffin",
"Monthly Interpolation")) +
theme_bw() +
theme(text = element_text(family = "Segoe UI Semilight",
size = text_size),
plot.title = element_text(hjust = 0.5),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "top")
- Monthly evaporation stable isotope measurements must be estimated since
it cannot be directly measured. We estimate $\delta^{18}O_{E}$ as follows
(Krabbenhoft et al., 1990, Eq. 5):
\begin{align}
\tag{5}
\delta^{18}O_{E} = \frac{(1/\alpha) * \delta^{18}O_{L} - h*\delta^{18}O_{A} -
\epsilon}{1 - h +10^{-3}\Delta \epsilon}
\end{align}
The relative humidity normalized to the temperature of the surface water (-),
$h$, is calculated as follows (Mook, 2000, Eq. 1.8):
\begin{align}
\tag{6}
h = \frac{RH}{100}*\frac{es_{A}}{es_{L}}
\end{align}
where $RH$ is the relative humidity (%), $es_{A}$ is the saturated vapor
pressure for the air (kPa), and $es_{L}$ is the saturated vapor pressure for
the lake (kPa). Saturated vapor pressure is calculated based on Allen et al. (1998)
using air temperature (degrees C) or lake temperature (degrees C) for $T$ as
appropriate:
\begin{align}
\tag{7}
es = 0.6108*exp(\frac{17.27T}{237 + T})
\end{align}
The equilibrium isotope fractionation factor (-), $\alpha$, is calculated
based on Gibson et al. (2016), where $T$ is the lake surface temperature
(K). Note that this formulation represents the fraction of isotopes in liquid
vs. the fraction in vapor (i.e., L/V form, $\alpha > 1$). Others, including
Krabbenhoft et al. (1990), use the inverse form of $\alpha$ (i.e., V/L form,
$\alpha < 1$) in their equations. In all equations presented in this
document, we consistently use the L/V form for alpha. If you calculate the
V/L form of alpha, simply replace $\alpha$ with $\frac{1}{\alpha}$ in
the equations presented in this document.
\begin{align}
\tag{8}
\alpha &= exp(-7.685*10^{-3} + \frac{6.7123}{T} - \frac{1666.4}{T^{2}} + \frac{350410}{T^{3}})
\end{align}
The kinetic fractionation factor (-), $\Delta \epsilon$, is calculated based
on Krabbenhoft et al. (1990):
\begin{align}
\tag{9}
\Delta \epsilon = 14.3*(1 - h)
\end{align}
The total fractionation factor (-), $\epsilon$, combines both the isotopic
fractionation factor and the kinetic fractionation factor (Krabbenhoft et
al., 1990):
\begin{align}
\tag{10}
\epsilon = 1000(1 - \frac{1}{\alpha}) + \Delta \epsilon
\end{align}
Lastly, the isotopic composition of the atmosphere, $\delta^{18}O_{A}$, is
calculated based on Gibson et al. (2016):
\begin{align}
\tag{11}
\delta^{18}O_{A} &= \frac{\delta^{18}O_{P} - k\epsilon^{+}}{1 + 10^{-3}k*\epsilon^{+}}
\end{align}
where $k$ can vary from 0.5 for highly seasonal climates to 1.0 for
non-seasonal climates but is set to 1.0 here and $\epsilon^{+}$ is defined
as:
\begin{align}
\tag{12}
\epsilon^{+} &= (\alpha - 1)*1000
\end{align}
plot_input(df = inputs,
ydata = "d18O_evap",
ylabel = expression(delta^18*'O'[evap]),
yscales = scales::number,
geom_type = "point",
text_size = text_size)
With the above information, we can then calculate groundwater inflow using Eq. 3.
plot_input(df = outputs,
ydata = "GWin_m3",
ylabel = expression(Groundwater~Inflow~(m^{3})),
inputs_only = FALSE,
text_size = text_size)
Groundwater Outflow {#GWout}
Finally, we calculate groundwater outflow as the only remaining term of the
water balance per Eq. 4.
plot_input(df = outputs,
ydata = "GWout_m3",
ylabel = expression(Groundwater~Outflow~(m^{3})),
inputs_only = FALSE,
text_size = text_size)
References
Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. (1998). Crop
evapotranspiration: Guidelines for computing crop water requirements. Rome: FAO.
Retrieved from http://www.fao.org/docrep/X0490E/x0490e00.htm.
Gibson, J.J., S.J. Birks, and Y. Yi. 2016. Stable isotope mass balance of lakes:
a contemporary perspective. Quaternary Science Reviews, 131:316-328.
https://doi.org/10.1016/j.quascirev.2015.04.013
Krabbenhoft, D. P., C. J. Bowser, M. P. Anderson, and J. W. Valley. (1990).
Estimating Groundwater Exchange with Lakes: 1. The Stable Isotope Mass Balance
Method. Water Resources Research, 26(10):2445-2453.
https://doi.org/10.1029/WR026i010p02445
McJannet, D. L., Webster, I. T., Stenson, M. P., and Sherman, B.S. (2008).
Estimating open water evaporation for the Murray-Darling Basin. A report to the
Australian Government from the CSIRO Murray-Darling Basin Sustainable Yields
Project, CSIRO, Australia, 50 pp. Retrieved from
http://www.clw.csiro.au/publications/waterforahealthycountry/mdbsy/technical/U-OpenWaterEvaporation.pdf.
McMahon, T. A., Peel, M. C., Lowe, L., Srikanthan, R., and McVicar, T. R.
(2013). Estimating actual, potential, reference crop and pan evaporation using
standard meteorological data: a pragmatic synthesis, Hydrol. Earth Syst. Sci.,
17, 1331–1363. https://doi.org/10.5194/hess-17-1331-2013.
Mook, W.G. (ed.) 2000. Environmental Isotopes in the Hydrologic Cycle: Volume
III: Surface Water. UNESCO. Paris, France.
Session Info
The version of R and the versions of packages in use when this R Markdown file
was created are listed here:
devtools::session_info()
cvoter/isoH2Obudget documentation built on March 29, 2020, 11:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Libraries library(devtools) library(CSLSiso) library(CSLSdata) library(reshape2) # for melting library(ggplot2) # for plotting library(cowplot) # for subplots library(extrafont) # for fonts library(dplyr) # for %>% pipes library(lubridate) library(knitr) library(stringr) library(scales) # for :: # Parameters text_size <- 12 lakes <- c("Pleasant", "Long", "Plainfield") # Speed up processing by getting combinded inputs and outputs once here. inputs <- NULL outputs <- NULL for (lake in lakes){ # Inputs these_inputs <- runall_csls_inputs(lake) these_inputs$lake <- lake inputs <- rbind(inputs, these_inputs) # Outputs these_outputs <- runall_csls_budget(lake) these_outputs$lake <- lake outputs <- rbind(outputs, these_outputs) } inputs$lake <- factor(inputs$lake, levels = lakes) outputs$lake <- factor(outputs$lake, levels = lakes)
Overview {#top}
This water budget model uses 1) observed changes in lake levels, 2) observed precipitation, 3) calculated lake evaporation, and 4) stable isotope measurements of precipitation, evaporation, lake water, and inflowing groundwater to calculate monthly groundwater flow into and out of a lake.
The lakes in the Central Sands Lakes Study (Pleasant Lake, Long Lake, and Plainfield Lake) are all seepage lakes, a type of lake that does not have inflowing streams, outflowing streams, or appreciable surface runoff. The water budget for a seepage lake is as follows:
\begin{align} \tag{1} \frac{dV}{dt} = P - E + GW_{in} - GW_{out} \end{align}
where $\frac{dV}{dt}$ is the change in lake volume (m^3^), $P$ is the precipitation (m^3^), $E$ is the lake evaporation (m^3^), $GW_{in}$ is the groundwater flowing in to the lake (m^3^), and $GW_{out}$ is the groundwater flowing out of the lake.
Change in Lake Volume is directly measured using a water level logger at each lake. Daily lake levels elevations are converted into daily lake volumes using lake elevation-volume relationships developed from lake bathymetry maps.
Precipitation is directly measured at a nearby weather station. It is converted to a volume by multiplying the total precipitation (mm) by the average lake area during the month of interest. Average lake area is determined using daily lake level elevation from water level loggers and the lake elevation-area relationships developed from bathymetry maps.
Lake Evaporation is calculated based on the lake energy-budget approach outlined in McJannet et al. (2008) as described in McMahon et al. (2013). This approach uses:
- Measured incoming solar radiation, wind speed, relative humidity, and air temperature from a nearby weather station
- Lake surface temperature from continuous temperature loggers
- Area and depth of the lake from daily lake stage measurements and elevation-area relationships
Evaporation as a depth (Eq. 2) is summarized at a monthly timestep and converted to a volume using the average lake area during that month, as in the precipitation calculations.
\begin{align} \tag{2} E = \frac{\Delta_{w}(R_{n} - G_{w}) + \frac{606024\rho_{a}c_{a}(es_{w} - es_{a})}{r_{a}}}{\lambda*(\Delta_{w} + \gamma)} \end{align}
where $E$ is lake evaporation (mm/day), $\Delta_{w}$ is the slope of the saturation water vapour curve at water temperature (kPa/deg C), $R_{n}$ is the net radiation (MJ/m^2^/day), $G_{w}$ is the change in heat storage in the water body (MJ/m^2/day), $\rho_{a}$ is the density of air (kg/m^3^), $c_{a}$ is the specific heat of air (MJ/kg/K), $es_{w}$ is the saturation vapour pressure at water temperature (kPa), $es_{a}$ is the saturation vapour pressure of the air (kPa), $r_{a}$ is the aerodynamic resistence (s/m), $\lambda$ is the latent heat of vapourisation (MJ/kg), and $\gamma$ is the psychormetric constant (kPa/deg C).
Groundwater inflow is calculated using monthly changes in lake volume, precipitation, lake evaporation, and stable isotope measurements based on the approach in Krabbenhoft et al. (1990) and modified based on Mook (2000) to account for non-steady state conditions in lake levels:
\begin{align} \tag{3} GW_{in} = \frac{P(\delta^{18}O_{L} - \delta^{18}O_{P}) + E(\delta^{18}O_{E} - \delta^{18}O_{L}) + V*d\delta^{18}O_{L}/dt}{\delta^{18}O_{GW_{in}} - \delta^{18}O_{L}} \end{align}
where $GW_{in}$ is the groundwater inflow (m^3^), $P$ is precipitation (m^3^), $E$ is evaporation (m^3^), V is the lake volume (m^3^) and $\delta^{18}O_{x}$ is the relative isotopic composition in units of per mil relative to a known standard. $\delta^{18}O_{L}$ is measured at the lake, $\delta^{18}O_{P}$ is measured from precipitation, $\delta^{18}O_{E}$ is estimated for evaporation, and $\delta^{18}O_{GW_{in}}$ is measured at inflowing groundwater wells.
Groundwater outflow can then be calculated as the only remaining term in the water balance:
\begin{align} \tag{4} GW_{out} = GW_{in} + P - E - \frac{dV}{dt} \end{align}
Detailed Approach
Input data for the Central Sands Lake Study sites (Pleasant Lake, Long Lake, and
Plainfield Lake) are saved in the CSLSdata R data package. The CSLSdata
package includes a data-raw directory where raw data is stored together with
the scripts used to clean raw data and save them as Rda files. Detailed
information on the source of each dataset is captured in the CSLSdata manual.
This is accessible as a pdf or via RStudio (help(package = CSLSdata)).
Below, we describe all additional assumptions and calculations used to transform
CSLSdata Rda files from their form in the CSLSdata package into the monthly
inputs required for monthly water balance calculations.
Change in Lake Volume {#dV}
We obtain daily lake elevation data from water level loggers at USGS gages at all three lakes in the Central Sands Lakes Study.
plot_lake_level(text_size = text_size)
We also have the elevation-volume relationships for each lake from lake bathymetry data (calculated in ArcGIS).
plot_elev_area_vol(ydata = "vol_m3", ylabel = expression(Lake~Volume~(m^{3})), text_size = text_size)
We convert these daily lake level elevations to daily lake volumes using these elevation-volume relationships.
plot_lake_level(ydata = "vol_m3", ylabel = expression(Lake~Volume~(m^{3})), yscales = scales::scientific, text_size = text_size)
The change in lake volume for a given month is then the difference between the lake volume on the last day of the month and the lake volume on the first day of the month.
plot_input(df = inputs, ydata = "dV_m3", ylabel = expression(Change~Lake~Volume~(m^{3})), geom_type = "col", text_size = text_size)
Precipitation {#precip}
We obtain hourly precipitation depth (mm) from the Hancock Agricultural Research Station (station id: hck). The station is located in Hancock, WI (location: 44.1188, -89.533, elevation: 241m), approximately 8 miles from Plainfield Lake, 8.5 miles from Long Lake, and 14.5 miles from Pleasant Lake.
weather <- CSLSdata::weather monthly_weather <- weather %>% group_by(date = floor_date(.data$date, unit = "month")) %>% summarise(P = sum(.data$P)) ggplot(monthly_weather, aes(x = date, y = P)) + geom_col() + scale_x_datetime(date_breaks = "3 months", date_labels = "%b '%y") + scale_y_continuous(labels = scales::scientific) + labs(x = "", y = "Monthly Precipitation (mm)") + theme_bw() + theme(text = element_text(family = "Segoe UI Semilight", size = text_size), plot.title = element_text(hjust = 0.5), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
To convert precipitation as a depth to precipitation as a volume, We use the daily lake elevations at the USGS gages at each lake.
plot_lake_level(text_size = text_size)
We also use the elevation-area relationships for each lake from lake bathymetry data (calculated in ArcGIS).
plot_elev_area_vol(ydata = "area_m2", ylabel = expression(Lake~Area~(m^{2})), text_size = text_size)
We convert daily lake level elevations to daily lake areas using these elevation-area relationships.
plot_lake_level(ydata = "area_m2", ylabel = expression(Lake~Area~(m^{2})), yscales = scales::scientific, text_size = text_size)
Then, we calculate the mean monthly lake area.
plot_input(df = inputs, ydata = "mean_area_m2", ylabel = expression(Lake~Area~(m^{2})), geom_type = "point", facet_scales = "free_y", text_size = text_size)
By multiplying monthly precipitation as a depth (mm) by the mean monthly lake area (m^2^), we obtain monthly precipitation as a volume (m^3^) for each lake.
plot_input(df = inputs, ydata = "P_m3", ylabel = expression(Precipitation~(m^{3})), geom_type = "col", text_size = text_size)
Lake Evaporation {#evap}
All lake evaporation calculations are contained in the CSLSevap package; more
details on this energy-balance approach can be obtained by searching the
documentation for that package. Broadly, the energy-budget lake evaporation
equations require information on:
- Lake location (latitude, longitude, elevation, latitude of time zone)
- Lake size (area, depth)
- Lake surface temperature (daily values, taken from surface-most continuous temperature logger)
- Meteorological parameters (daily relative humidity, air temperature, incoming solar radiation, wind speed from nearby weather station)
- Heat absorption and reflection parameters (e.g., albedo of water, density and specific heat of air)
These lake evaporation equations yeild daily lake evaporation (mm/day), which is converted to monthly lake evaporation volume by summing over each month and multiplying by the average monthly lake area, as in the precipitation calculations.
plot_input(df = inputs, ydata = "E_m3", ylabel = expression(Evaporation~(m^{3})), geom_type = "col", text_size = text_size)
Groundwater Inflow {#GWin}
As indicated in Eq. 3, calculation of groundwater inflow requires:
-
Monthly precipitation, evaporation, and mean lake volume. These are calculated per the previous sections.
-
Monthly groundwater stable isotope measurements are derived from regular CSLS sampling at four groundwater wells near each lake. These measurements were taken approximately every other month. Upgradient wells representing inflowing groundwater are defined as wells where the mean daily difference between groundwater levels and lake levels in the 30-day period prior to sampling is greater than 1 cm (i.e., the precision of the water level logger). All upgradient well measurements from the same sampling event for a lake are averaged, then linearly interpolated to estimate values on the last day of each month.
-
Monthly lake stable isotope measurements are derived from regular CSLS sampling at each lake. These measurements were taken less frequently than groundwater measurements; approximately every other month during the ice-free period. All lake measurements are linearly interpolated to estimate values on the last day of each month during the ice-free period.
-
Monthly precipitation stable isotope measurements are derived from regular CSLS sampling at two precipitation collectors located at Hancock Agricultural Research Station. These measurements were taken approximately every month during the snow-free season. They are supplemented by monthly measurements taken by Maribeth Kniffin at the same location from May 2016 to April 2017. Since 1) we do not have access to the exact days that Maribeth Kniffin collected her measurements and 2) precipitation samples are an aggregate sample of all precipitation collected between sample dates (i.e., they are not a snapshot in time), we do not perform any linear interpolation on these samples. Instead, we assume all precipitation measurements represent precipitation on the last day of the month in which they were collected. Maribeth Kniffin's data is only used for months when we do not have CSLS-collected data.
lakes <- c("Pleasant", "Long", "Plainfield") monthly_isotopes <- NULL for (lake in lakes){ isotopes <- CSLSdata::isotopes[[lake]] dictionary <- CSLSdata::dictionary[[lake]] lake_levels <- CSLSdata::lake_levels[[lake]] gw_levels <- CSLSdata::gw_levels[[lake]] # Measured values tmp1 <- iso_site_type(isotopes, dictionary, lake_levels, gw_levels) tmp1 <- tmp1 %>% group_by(floor_date(.data$date, unit = "month"), variable = .data$site_type) %>% filter(is.na(.data$site_type) == FALSE, .data$site_type != "", .data$site_type != "downgradient", .data$site_type != "deep") %>% summarise(date = floor_date(mean(.data$date), unit = "day"), value = mean(.data$d18O)) %>% ungroup() %>% select(.data$date, .data$variable, .data$value) tmp1$variable <- tmp1$variable %>% str_replace("upgradient","Groundwater") %>% str_replace("lake","Lake") %>% str_replace("precipitation","Precipitation") tmp1$lake <- lake tmp1$measurement <- "actual" # monthly values tmp2 <- summarise_isotopes(isotopes, dictionary, find_timeseries(isotopes), lake_levels, gw_levels) %>% select(date, d18O_GWin, d18O_lake, d18O_pcpn) %>% melt(id.vars = "date") tmp2$variable <- tmp2$variable %>% str_replace("d18O_GWin","Groundwater") %>% str_replace("d18O_lake","Lake") %>% str_replace("d18O_pcpn","Precipitation") tmp2$lake <- lake tmp2$measurement <- "monthly" tmp3 <- tmp1 %>% filter(variable == "Precipitation") tmp3 <- fill_timeseries_gaps(tmp3, find_timeseries(isotopes)) kniffin_isotopes <- CSLSdata::kniffin$isotopes for (i in 1:nrow(tmp3)) { this_month <- month(tmp3$date[i]) this_i <- which(month(tmp3$date) == this_month) tmp3$value[i] <- mean(tmp3$value[this_i], na.rm = TRUE) tmp3$measurement[i] <- "actual" tmp3$value[is.nan(tmp3$value)] <- NA this_kniffin <- kniffin_isotopes %>% filter(month(.data$date) == this_month, .data$site_id == "PRECIP") if (nrow(this_kniffin) == 0) { kniffin_d18O <- NA } else { kniffin_d18O <- mean(this_kniffin$d18O, na.rm = TRUE) } if (is.na(tmp3$value[i])){ tmp3$value[i] <- kniffin_d18O tmp3$measurement[i] <- "kniffin" } } tmp3$lake <- lake tmp3$variable <- "Precipitation" monthly_isotopes <- rbind(monthly_isotopes, tmp1, tmp2, tmp3) } monthly_isotopes$lake <- factor(monthly_isotopes$lake, levels = lakes) ggplot(monthly_isotopes %>% arrange(measurement), aes(x = date, y = value, color = measurement, shape = measurement)) + geom_point() + scale_x_datetime(date_breaks = "4 months", date_labels = "%b '%y") + facet_grid(variable~lake) + labs(x = "", y = expression(delta^18*'O'), title = "") + scale_color_manual(name = "", values = c("red", "purple", "grey70"), labels = c("Measured by CSLS", "Measured by M. Kniffin", "Monthly Interpolation")) + scale_shape_manual(name = "", values = c(19, 15, 17), labels = c("Measured by CSLS", "Measured by M. Kniffin", "Monthly Interpolation")) + theme_bw() + theme(text = element_text(family = "Segoe UI Semilight", size = text_size), plot.title = element_text(hjust = 0.5), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "top")
- Monthly evaporation stable isotope measurements must be estimated since it cannot be directly measured. We estimate $\delta^{18}O_{E}$ as follows (Krabbenhoft et al., 1990, Eq. 5):
\begin{align} \tag{5} \delta^{18}O_{E} = \frac{(1/\alpha) * \delta^{18}O_{L} - h*\delta^{18}O_{A} - \epsilon}{1 - h +10^{-3}\Delta \epsilon} \end{align}
The relative humidity normalized to the temperature of the surface water (-), $h$, is calculated as follows (Mook, 2000, Eq. 1.8): \begin{align} \tag{6} h = \frac{RH}{100}*\frac{es_{A}}{es_{L}} \end{align}
where $RH$ is the relative humidity (%), $es_{A}$ is the saturated vapor pressure for the air (kPa), and $es_{L}$ is the saturated vapor pressure for the lake (kPa). Saturated vapor pressure is calculated based on Allen et al. (1998) using air temperature (degrees C) or lake temperature (degrees C) for $T$ as appropriate: \begin{align} \tag{7} es = 0.6108*exp(\frac{17.27T}{237 + T}) \end{align}
The equilibrium isotope fractionation factor (-), $\alpha$, is calculated based on Gibson et al. (2016), where $T$ is the lake surface temperature (K). Note that this formulation represents the fraction of isotopes in liquid vs. the fraction in vapor (i.e., L/V form, $\alpha > 1$). Others, including Krabbenhoft et al. (1990), use the inverse form of $\alpha$ (i.e., V/L form, $\alpha < 1$) in their equations. In all equations presented in this document, we consistently use the L/V form for alpha. If you calculate the V/L form of alpha, simply replace $\alpha$ with $\frac{1}{\alpha}$ in the equations presented in this document.
\begin{align} \tag{8} \alpha &= exp(-7.685*10^{-3} + \frac{6.7123}{T} - \frac{1666.4}{T^{2}} + \frac{350410}{T^{3}}) \end{align}
The kinetic fractionation factor (-), $\Delta \epsilon$, is calculated based on Krabbenhoft et al. (1990): \begin{align} \tag{9} \Delta \epsilon = 14.3*(1 - h) \end{align}
The total fractionation factor (-), $\epsilon$, combines both the isotopic fractionation factor and the kinetic fractionation factor (Krabbenhoft et al., 1990): \begin{align} \tag{10} \epsilon = 1000(1 - \frac{1}{\alpha}) + \Delta \epsilon \end{align}
Lastly, the isotopic composition of the atmosphere, $\delta^{18}O_{A}$, is calculated based on Gibson et al. (2016): \begin{align} \tag{11} \delta^{18}O_{A} &= \frac{\delta^{18}O_{P} - k\epsilon^{+}}{1 + 10^{-3}k*\epsilon^{+}} \end{align}
where $k$ can vary from 0.5 for highly seasonal climates to 1.0 for non-seasonal climates but is set to 1.0 here and $\epsilon^{+}$ is defined as: \begin{align} \tag{12} \epsilon^{+} &= (\alpha - 1)*1000 \end{align}
plot_input(df = inputs, ydata = "d18O_evap", ylabel = expression(delta^18*'O'[evap]), yscales = scales::number, geom_type = "point", text_size = text_size)
With the above information, we can then calculate groundwater inflow using Eq. 3.
plot_input(df = outputs, ydata = "GWin_m3", ylabel = expression(Groundwater~Inflow~(m^{3})), inputs_only = FALSE, text_size = text_size)
Groundwater Outflow {#GWout}
Finally, we calculate groundwater outflow as the only remaining term of the water balance per Eq. 4.
plot_input(df = outputs, ydata = "GWout_m3", ylabel = expression(Groundwater~Outflow~(m^{3})), inputs_only = FALSE, text_size = text_size)
References
Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. (1998). Crop evapotranspiration: Guidelines for computing crop water requirements. Rome: FAO. Retrieved from http://www.fao.org/docrep/X0490E/x0490e00.htm.
Gibson, J.J., S.J. Birks, and Y. Yi. 2016. Stable isotope mass balance of lakes: a contemporary perspective. Quaternary Science Reviews, 131:316-328. https://doi.org/10.1016/j.quascirev.2015.04.013
Krabbenhoft, D. P., C. J. Bowser, M. P. Anderson, and J. W. Valley. (1990). Estimating Groundwater Exchange with Lakes: 1. The Stable Isotope Mass Balance Method. Water Resources Research, 26(10):2445-2453. https://doi.org/10.1029/WR026i010p02445
McJannet, D. L., Webster, I. T., Stenson, M. P., and Sherman, B.S. (2008). Estimating open water evaporation for the Murray-Darling Basin. A report to the Australian Government from the CSIRO Murray-Darling Basin Sustainable Yields Project, CSIRO, Australia, 50 pp. Retrieved from http://www.clw.csiro.au/publications/waterforahealthycountry/mdbsy/technical/U-OpenWaterEvaporation.pdf.
McMahon, T. A., Peel, M. C., Lowe, L., Srikanthan, R., and McVicar, T. R. (2013). Estimating actual, potential, reference crop and pan evaporation using standard meteorological data: a pragmatic synthesis, Hydrol. Earth Syst. Sci., 17, 1331–1363. https://doi.org/10.5194/hess-17-1331-2013.
Mook, W.G. (ed.) 2000. Environmental Isotopes in the Hydrologic Cycle: Volume III: Surface Water. UNESCO. Paris, France.
Session Info
The version of R and the versions of packages in use when this R Markdown file was created are listed here:
devtools::session_info()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.