# In cvoter/isoH2Obudget: Lake Water Budget Using Stable Isotopes

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 ::

text_size = 12


## 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 roxygen documentation files for this package. For example, to access information on the lake level dataset, load the package into RStudio using library(CSLSdata) and try ?lake_levels.

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.

lake_level_plot <- function(ydata, ylabel, title = "",
lakes = c("Pleasant", "Long", "Plainfield")){
lake_levels <- NULL
for (lake in lakes){
tmp           <- CSLSdata::lake_levels[[lake]]
tmp           <- tmp[c("date", ydata)]
colnames(tmp) <- c("x", "y")
tmp$lake <- lake lake_levels <- rbind(lake_levels, tmp) } lake_levels$lake <- factor(lake_levels$lake, levels = lakes) plot_obj <- ggplot(lake_levels) + geom_line(aes(x = x, y = y)) + scale_x_datetime(date_breaks = "6 months", date_labels = "%b '%y") + labs(x = "", y = ylabel, title = title) + facet_wrap(~lake, scales = "free_y", ncol = 3) + 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()) if (ydata != "level_m") { plot_obj <- plot_obj + scale_y_continuous(labels = scales::scientific) } return(plot_obj) } elev_area_vol_plot <- function(ydata, ylabel, title = "", lakes = c("Pleasant", "Long", "Plainfield")){ elev_area_vol <- NULL for (lake in lakes){ tmp <- CSLSdata::elev_area_vol[[lake]] tmp <- tmp[c("elev_m", ydata)] colnames(tmp) <- c("x", "y") tmp$lake      <- lake
elev_area_vol <- rbind(elev_area_vol, tmp)
}
elev_area_vol$lake <- factor(elev_area_vol$lake, levels = lakes)

plot_obj <- ggplot(elev_area_vol) +
geom_line(aes(x = x, y = y)) +
facet_wrap(~lake, scales = "free") +
scale_y_continuous(labels = scales::scientific) +
labs(x = "Elevation (m)", y = ylabel, title = title) +
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())

return(plot_obj)
}

inputs_plot <- function(ydata, ylabel, geom_type = "col", title = "",
lakes = c("Pleasant", "Long", "Plainfield")){
inputs <- NULL
for (lake in lakes){
tmp           <- runall_csls_inputs(lake)
tmp           <- tmp[c("date", ydata)]
colnames(tmp) <- c("x", "y")
tmp$lake <- lake inputs <- rbind(inputs, tmp) } inputs$lake <- factor(inputs$lake, levels = lakes) # Specify plot type (column or point) if (geom_type == "col") { plot_obj <- ggplot(inputs, aes(x = x, y = y)) + geom_col() } else if (geom_type == "point") { plot_obj <- ggplot(inputs, aes(x = x, y = y)) + geom_point() } plot_obj <- plot_obj + scale_x_datetime(date_breaks = "6 months", date_labels = "%b '%y") + facet_wrap(~lake) + scale_y_continuous(labels = scales::scientific) + labs(x = "", y = ylabel, title = title) + 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()) return(plot_obj) } iso_plot <- function(ydata = "d18O", ylabel = expression(delta^18*'O'), title = "", 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 != "") %>% 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("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"

monthly_isotopes <- rbind(monthly_isotopes, tmp1, tmp2)
}

plot_obj <- ggplot(monthly_isotopes,
aes(x = date, y = value, color = measurement)) +
geom_point() +
scale_x_datetime(date_breaks = "6 months",
date_labels = "%b '%y") +
facet_wrap(lake~variable) +
labs(x = "", y = ylabel, title = title) +
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())

return(plot_obj)
}


### 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.

lake_level_plot("level_m", "Lake Elevation (m)")


We also have the elevation-volume relationships for each lake from lake bathymetry data (calculated in ArcGIS).

elev_area_vol_plot("vol_m3", expression(Lake~Volume~(m^{3})))


We convert these daily lake level elevations to daily lake volumes using these elevation-volume relationships.

lake_level_plot("vol_m3", expression(Lake~Volume~(m^{3})))


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.

inputs_plot("dV_m3", expression(Change~Lake~Volume~(m^{3})))


### 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.

lake_level_plot("level_m", "Lake Elevation (m)")


We also use the elevation-area relationships for each lake from lake bathymetry data (calculated in ArcGIS).

elev_area_vol_plot("area_m2", expression(Lake~Area~(m^{2})))


We convert daily lake level elevations to daily lake areas using these elevation-area relationships.

lake_level_plot("area_m2", expression(Lake~Area~(m^{2})))


Then, we calculate the mean monthly lake area.

inputs_plot("mean_area_m2", expression(Lake~Area~(m^{2})), "point")


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.

inputs_plot("P_m3", expression(Precipitation~(m^{3})), "col")


### Lake Evaporation {#evap}

All lake evaporation calculations are contained in the H2Oevap 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.

inputs_plot("E_m3", expression(Evaporation~(m^{3})), "col")


### Groundwater Inflow {#GWin}

As indicated in Eq. 3, calculation of groundwater inflow requires:

1. Monthly precipitation, evaporation, and mean lake volume. These are calculated per the previous sections.

2. 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.

3. 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.

4. 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") %>% 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("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 <- CSLSiso::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 kniffin_month <- kniffin_isotopes %>% filter(month(.data$date) == this_month)
if (is.na(tmp3$value[i])){ tmp3$value[i]       <- kniffin_month$d18O_pcpn 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")

1. 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{\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 Ozaydin et al. (2001), where $t$ is the lake surface temperature (K):

\begin{align} \tag{8} \alpha &= e^{\frac{-2.0667 - \frac{415.6}{t} + (\frac{1137}{t^2})*10^3}{1000}} \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}

### Groundwater Outflow {#GWout}

devtools::session_info()