Black_Ash_Wetland_Response_to_Future_Climate_Conditions.Rmd

---
author: Joe Shannon
title: Black Ash Wetland Response to Future Climate Conditions
date: "`r format(Sys.Date(), '%B %d, %Y')`"
output:
  bookdown::word_document2:
    reference_docx : output/dissertation_chapter/climate_change_impacts.docx
    fig_caption: TRUE
    number_sections: FALSE
    toc: FALSE
header-includes:
- \usepackage{float}
bibliography: resources/references.bib
---

```{r setup, echo = FALSE, message=FALSE, warning=FALSE}

knitr::opts_chunk$set(echo = FALSE, cache = FALSE)

```

Abstract
Future climate scenarios will likely lead to reduced connectivity between black
ash headwater wetlands and stream networks. 


## Black Ash Wetland Response to Future Climate Conditions

### Introduction
Shifts in temperature, and precipitation (timing, frequency, and quantity)
associated with climate change are already impacting global wetland ecology and
will intensify into the future [@burkett-2000;
@moomaw-2018]. In general, the Great Lakes region of North America (Figure
\@ref(fig:study-map)) is expect to see declining summer precipitation while total
annual precipitation stays stable or increases [@byun-2018;
@hayhoe-2010]. These projected changes result in increased precipitation in the
spring and winter months. Winter precipitation will also experience a phase
change, seeing a reduction in the proportion of precipitation as snow and
increased occurrence of rain-on-snow melt events in the winter and spring
[@notaro-2015]. Locally conditions will vary and our study area (see below) is 
projected in numerous downscaling scenarios across multiple GCMs to see a smaller 
than regional decline or even slight increases in summer precipitation [@pierce-2014, 
@hegewisch-2021, @byun-2018 Figure 7]. Taken 
together the shift in precipitation timing, the reduction of total snowmelt, and 
an earlier snowmelt will reshape the annual hydrologic budget. 
These changes will alter wetland hydroperiods, with water availability
increasing in the winter and spring and decreasing in the summer and early fall.
The impact of reduced or steady precipitation in the summer will be compounded by
increased evaporative demand [@byun-2018]. Summer temperatures are expected to
increase 2-8C° [@byun-2018; @hayhoe-2010], and without a commensurate 
increase in summer precipitation will lead to more frequent drought conditions 
and increased water stress on ecosystems.

Black ash (*Fraxinus nigra* Marsh) is an important hardwood component of many
forested wetlands in the northern United States and southern Canada. Wetlands 
with a large or majority black ash component face the same climate change 
impacts as other wetlands in the region. They also face loss of the 
major canopy species due to an invasive insect, Emerald Ash Borer (*Agrilus
planipennis* Fairmaire, EAB). Emerald ash borer was first found in the United
States in southeast Michigan in 2002 [@haack-2002]. It is known to infest and
cause high mortality in all ash native to North America [@herms-2014]. As of this
writing it is present in 35 US states and 5 Canadian provinces.
High mortality and the importance of ash in regional forested wetlands has led
to research on the impacts of their loss and strategies to mitigate those
impacts. @iverson-2016 evaluated the potential replacement species for black ash 
in the context of habitat availability, species migration, and replacement
species susceptibility to climate change impact, providing a useful resource for
climate-informed species replacements for black ash. However, black ash grows in
a range of geomorphic settings and local site conditions can have a strong
influence on hydrology, forest structure, and plant community among wetlands
[@kolka-2018]. Considering both @iverson-2016 and @kolka-2018 we can identify
four necessary components to evaluate mitigation efforts in black ash wetlands:
1) climate-informed species selection, 2) species tolerance of local site
conditions, 3) site conditions following EAB infestation, and 4) site conditions
in a future climate.

Previous and ongoing work in Michigan and Minnesota has assessed mitigation in
light of three of the four components of the combined impact of climate change
and EAB on black ash wetlands. Researchers in Michigan [@bolton-2018],
Minnesota [@looney-2015], and Wisconsin [@bolton-2018], planted seedlings to
evaluate potential canopy species at the wetland rather than landscape level.
The plantings in @bolton-2018 and @looney-2015 took place under simulated EAB
infestation, where the seedlings were subjected to adverse conditions due to
the hydrologic impact of the loss of black ash [@slesak-2014;
@vangrinsven-2017], as well as the increased competition from herbaceous growth
under increased light conditions [@looney-2016; @looney-2017a; @davis-2017].

The researchers addressed 3 of the 4 components necessary for mitigation
efforts in black ash wetlands: 1) some planted species were at the northern
edge of their current range and evaluated 2) under present-day site conditions
and 3) under simulated-EAB site conditions. What is not possible in field
studies is addressing the interaction of EAB and climate change on site
conditions. Just as we have seen EAB impacts cascade through black ash
ecosystems affecting hydrology, plant communities, and nutrient cycling, we can
expect the future climate-driven changes to hydrology to result in similar
cascades. Focusing on the future hydrologic characteristics of these wetlands
will inform us on the impacts to other functions, as hydrology is the critical
control on wetland ecosystems [@brinson-1993]. In order to understand future
hydrologic conditions we have derive wetland water level models and evaluate
under future climate scenarios.

Simulated daily time-step weather conditions are available for various future
climate scenarios to serve as model inputs, providing realizations of 30-year
climate periods. However, previous studies have demonstrated high-interannual
variation in seasonal wetland water level drawdown and rebound in black ash
wetlands [@vangrinsven-2017]. Estimating the behavior of a highly variable
system with 30 years of daily climate projects could lead to biased or
high-variance results making it difficult to draw proper conclusions about the
combined impacts of climate change and EAB. A 30-year period provides a
representation of normal climate variation, but the weather sequences
underlying a certain climate can take many more forms than the observed (or
simulated) 30 annual weather sequences. Therefore, simulation studies of
alternative scenarios or future conditions require more weather sequences to
quantify the most likely response and the distribution of possible responses.
Stochastic weather generators (SWGs) provide a tool for generating synthetic
time series of weather that simulate conditions under an observed or projected
climatology [@wilks-1999]. SWGs have also been used as a downscaling technique
for global or regional climate models [@wilks-2012; @verdin-2015;
@verdin-2019]. One class of SWG known as Richardson [@wilks-2012] are built
from parametric representations fit to climate conditions using multiple
generalized linear models and varied statistical distributions. This approach
can be implemented using any observed (or simulated) weather series that define
particular climatic conditions. This characteristic makes them well suited to
compare multiple climate scenarios without the overhead cost of executing
complete GCM or regional climate model (RCM) runs.

We have performed simulation experiments combining observed wetland hydrology
and synthetic weather sequences to quantify the interactions of EAB and future
climate scenarios. We developed wetland hydrology models for current canopy
conditions (ash-dominated forested wetlands), non-forested conditions
(herbaceous and shrub/scrub wetlands), and potential future canopy conditions
(forested wetlands under current co-dominant species). We evaluated each class
of model under two potential future climate scenarios for the end of the 21st
century (2070-2099). The two future scenarios are defined by a less sensitive
GCM (generally projects smaller magnitude regional climate impacts) under the
moderate relative concentration pathway (RCP 4.5), and a more sensitive GCM
(projects larger magnitude regional climate impacts) under the relative
concentration pathway (RCP 8.5) [@vanvuuren-2011]. The two-scenario ‘bookend’
approach provides a range of potential future conditions as opposed to a
multi-model ensemble approach which masks some of the uncertainty in potential
future conditions [@swanston-2016]. For an unbiased comparison between current
and future conditions both GCMs are used to evaluate the wetland models under
historic (1980-2009) climate conditions.

We expect that the interaction of hydrologic impacts of EAB and climate change
will result in a tempering of the two individual impacts in this region. While
simulation of post-EAB conditions have lead to increased water levels and
reduced drawdown rates in growing season, future climate conditions in the
region will result in reduced water availability during that same period. These
two opposing drivers should result in some moderation to both impacts.

### Methods

#### Study Area & Data

Wetland water levels measured at fourteen black ash-dominated wetlands in the
western Upper Peninsula of Michigan, USA were used to develop and evaluate our
wetland hydrologic models (Fig. \@ref(fig:study-map)). The wetlands range in
size from 0.29 to 1.19 ha and 30--80% of the basal area consists of black ash
with histosol soils over a confining layer located at an average depth of 118.8
cm [@davis-2017; @vangrinsven-2017]. The region has average minimum and maximum
annual temperatures of -11.3 °C and 18.2 °C and an average annual
precipitation of 101 cm for the climate period of 1980-2009 at the Bergland Dam
(46°35'13"N, 89°32'51"W) meteorological station [@arguez-2012]. Wetland
water levels were continuously monitored in 2" inner-diameter driven wells from
2012-2020 and logged every fifteen minutes using Solinst Levellogger Junior
pressure transducers (Solinst, Ontario, CA), with more details available in
[@vangrinsven-2017]. Barometric compensation was performed using data from
Solinst Levellogger Junior pressure transducers deployed at a subset of study
wetlands. Compensated water levels were corrected for temperature differentials
as in Shannon, et al (In Review). Following an initial control period of two
growing seasons, two-thirds of the wetlands in the study were treated to 
simulate the impacts of an EAB infestation. At one-third of the sites all ash 
trees greater than 1" in diameter were girdled and at the other half of the 
treatment sites all ash trees greater than 1" in diameter were hand-felled and 
left on site.

Daily precipitation and daily minimum/maximum temperatures from existing
meteorological stations were used as inputs to the wetland water level models
described below. From these input drivers we derived daily solar radiation, PET,
precipitation as snowfall and rain, and snowmelt. Precipitation records were
retrieved from the National Centers for Environmental Information Hourly
Precipitation Dataset [@hpdn] using the stations USC00201088 (Bruce Crossing,
MI), USC00204328 (Kenton, MI), USC00475352 (Mercer Ranger Station, WI),
USC00206215 (Ontonagon, MI), USC00476398 (Park Falls, WI), USC00476518 (Phelps,
WI), USC00476939 (Rainbow Reservoir Lake, Tomahawk, WI), USC00477140 (Rice
Reservoir, Tomahawk, WI), and USC00208680 (Watersmeet Fish Hatchery) and summed
to daily values. Daily minimum and maximum temperatures were taken from the
Global Historical Climatology Network (GHCND) dataset [@ghcnd; @menne-2012] for
the same stations, which were collocated stations for HPD and GHCND. Where data were
missing, values were filled using inverse distance weighting with data retrieved
from the Mesowest [@horel-2002] network stations BPLM4 (Baraga Plains, MI),
KTNM4 (Kenton, MI), PIEM4 (Pelkie, MI), WKFM4 (Wakefield, MI), WMTM4
(Watersmeet, MI). Solar radiation data were also retrieved from the listed
Mesowest stations. The solar radiation and temperature data were used to fit the
Bristow-Campbell method coefficients to estimate solar radiation from latitude
and daily temperature range using the PIEM4 Mesowest site [@bristow-1984;
@bojanowski-2016]. Bristow-Campbell solar radiation was used to calculate
potential evapotranspiration (PET) via the modified Hargreaves-Samani equation
[@hargreaves-2003]. Precipitation was partitioned into snowfall, rain, and
snowmelt inputs using the CemaNeige snow accounting routine (SAR)
[@valery-2014b]. The CemaNeige SAR is temperature-index based and accounts for
accumulation, snowpack cold content, and snowmelt through a thermal state
weighting coefficient and a degree-day melt coefficient. The CemaNeige SAR is
implemented in the R package `airGR` [@coron-2017].

```{r study-map, fig.cap = "Map of study site and meteorological stations used for data retreival. Coordinates in meters, UTM Zone 16N."}

knitr::include_graphics("output/figures/climate_study_map.png")

```

#### Wetland Hydrology Models

The link between hydrology drivers (rainfall, snowmelt, PET) and daily water
level response has been shown to be non-linear and varying with stage
[@loheideii-2008; @white-1932]. In wetland systems the relationship between the
driver magnitude and response magnitude has been termed the ecosystem specific
yield (E~Sy~) [@mclaughlin-2014] (see **Section #### in Chapter #### on
Pressure Transducer Uncertainty for more information about the development of
E~Sy~ over time**). E~Sy~ has previously been empirically derived as the ratio
between precipitation inputs and water level rise ($\frac{P}{\Delta\;WL}$)
[@mclaughlin-2014]. The relationship between empirical E~Sy~ and water level
can then be modeled to provide a continuous estimate of E~Sy~. Models fit to
the E~Sy~\~Water Level relationship include exponential [@mclaughlin-2014;
@watras-2017], quadratic [@mclaughlin-2014], and step-wise regression
[@mclaughlin-2014]. Identifying *E~Sy~* using the rainfall-rise ratio can be
unworkable in the presence of confounding hydrologic variation such as surface
water connectivity and low-frequency seasonal changes in water availability
[@watras-2017; @zhu-2011]. Both of these factors were present in our study 
wetlands and no clear relationship existed to model *E~Sy~* with the above 
model forms, requiring an alternative approach.

As an alternative we fit an inverse analog to $\frac{P}{\Delta\;WL}$, deriving
*E~Sy~* from the ratio of cumulative water availability and water level. We
began by fitting a quadratic curve to the relationship between a year-to-date
water availability index from the beginning of the growing season to the point
of minimum wetland water level (Figure \@ref(fig:esy-function)), where

$$WA_{YTD}\;=\;P_{YTD} + M_{YTD} - PET_{YTD},$$  

$WA_{YTD}$, $P_{YTD}$, $M_{YTD}$, and $PET_{YTD}$ are year-to-date water
availability, rainfall, snow melt, and potential evapotranspiration,
respectively. Empirical *E~Sy~* was then considered as the first derivative of the
quadratic curve, which has the form $\frac{\Delta\;WL}{\Delta\;WA}$. This
approach resulted in an asymptotic relationship between *E~Sy~* and water level,
suggesting agreement with the exponential forms in
@mclaughlin-2014 and @watras-2017. Defining *E~Sy~* as
$\frac{\Delta\;WL}{\Delta\;WA}$ provided additional information to fit the
relationship between *E~Sy~* and water level as we were not limited to days with
rainfall. Models for *E~Sy~* were fit using the year with the greatest water level
drawdown for each wetland to capture the widest range of *E~Sy~* variation. Each
wetland had the same basic form of *E~Sy~* function. The functions were fit using
a single hierarchical model with each of the coefficients allowed to vary
independently within sites. Models were fit using the `brms` package in R
[@burkner-2017; @burkner-2018]. The asymptotic structure of the *E~Sy~* function
requires that a lower bound be placed on *E~Sy~* predictions to avoid values less
than or equal to zero. The minimum value of *E~Sy~* was allowed to vary by wetland
and was fit along with other model parameters as described below.

```{r esy-function, fig.cap="Derivation of ecosystem specific yield for a single study wetland. Panel A shows the quadratic relationship between wetland water level and water balance (total liquid inputs minus potential evapotranspiration) for the drawdown period. Panel B shows the relationship between the derivative of the fitted line from Panel A, representing an emprical estimate of E~Sy~ and wetland water level."}

knitr::include_graphics("output/figures/scatter_plot_and_line_empirical_esy.png")

```


Wetland water levels were simulated on a continuous basis using daily inputs.
For each daily step, wetland water level change was determined as a function of
the current water level, daily rainfall (*R*), daily *PET*, daily snowmelt (*M*), and
estimated streamflow (*Q*) (Equation EQAAA; Table \@ref(tab:parameter-table)). Maximum
wetland water levels were estimated from the the control period as the mode of the
wetland water level record. Streamflow was assumed to occur whenever wetland
water levels were at or above the maximum water level. This is similar to the
approach in @mclaughlin-2019 where wetland surface water connectivity was
determined from wetland water level records. The current water level was used to
estimate $\hat_{E}~Sy~$ which was then used as a multiplier for rainfall, *PET*,
and snowmelt components. In addition each driver (*R*, *PET*, *M*, *Q*) had a 
coefficient (fitting described below). Snowmelt and precipitation were also fit
with a first order autoregressive filter to simulate slow flow contributions to
wetland water levels. Within a single day only the larger of *PET* or *P* was used 
as driver of water level change, accounting for suppressed transpiration in from
wet leaves.

$$
\begin{align*}
& WL_{t+1} = WL_{t} + f_{E_{Sy}}(WL_{t-1})*(\hat{\beta}_DD + \hat{\beta}_MM) + \hat{\beta}_Q(WL_{t} - WL_{max})\\
& \hat{\beta}_DD\;=\left\{R \le PET;\;\;\hat{\beta}_{PET}PET\atop{R > PET; \;\; \hat{\beta}_RR}\;\;\;\;\;\;\;\;\;\;\right.\\
& f_{E_{Sy}}(WL_{t-1})\;=\;a - (a - b) * e^{(c * WL_{t-1})} \\
& M_t = M_t + \phi_M\; M_{t-1} \\
& R_t = R_t + \phi_R\;R_{t-1}
\end{align*}
$$

```{r parameter-table, echo = FALSE, warning = FALSE, message=FALSE}

library(flextable)
library(ftExtra)
library(officer)

tab <- data.table::fread("Parameter | Definition  
t | Daily time step
WL | Wetland water level (relative to ground surface)
WL~max~ | Maximum wetland water level (relative to ground surface)
M, $\\hat{\\beta}_M$ | Snowmelt and snowmelt coefficient
Q, $\\hat{\\beta}_Q$ | Streamflow and streamflow coefficient
R, $\\hat{\\beta}_R$ | Rainfall and rainfall coefficient
PET, $\\hat{\\beta}_PET$ | Potential evapotranspiration and potential evapotranspiration coefficient
$f_{E_{Sy}}$ | Estimate of ecosystem specific yield
$\\phi_M$ | Autoregressive coefficient for snowmelt
$\\phi_R$ | Autoregressive coefficient for rainfall
a, b, c | Fitted coefficients for asymptotic regression", sep = "|")

flextable(tab) %>% 
  colformat_md() %>% 
  set_caption("Parameter notation and description for wetland water level models.",
              autonum = run_autonum(seq_id = "tab", bkm = "parameter-table")) %>% 
  autofit()

```

Training data for wetland water level models was selected for each wetland as
the year with the greatest water level drawdown within the control period. The
remaining control-period years were used for wetland model evaluation. Wetland
parameters were estimated using the L-BFGS-B algorithm via the R `optim()`
function [@rcoreteam-2019] minimizing the weighted root-mean square error of
the predictions. Weights were applied asymmetrically where days with observed
water levels greater than or equal to *WL<sub>max</sub>* were assigned the
original equal weight (1/n). As observed water levels decreased below
*WL<sub>max</sub>* weights increased as the square of the difference between
observed water level and *WL<sub>max</sub>*. To model impacted wetlands that
have not recovered to a fully forested state the observed treatment data from
the same site were used to reparameterize the control condition model.
Treatment-period training data were also selected to maximize water level
drawdown in the training set. The wetland models were reparameterized for only
$\hat{\beta}_{PET}$, $\hat{\beta}_{R}$ terms. Finally we simulated reforested
black ash wetlands under one set of potential future forest composition. For
the future forest conditions we assumed a mix of the current co-dominant
species, red maple and yellow birch, would become established with similar
stand basal area. Future forested simulations were performed using the same
parameters as the control period with a reduction of $\hat{\beta}_{PET}$. The
reduction in $\hat{\beta}_{PET}$ was a function of water level and current
proportion of site basal area as black ash:

$$\hat{\beta}_{PET-Adj} = \hat{\beta}_{PET}*[(1 - BA_{ash}) + BA_{ash} * [1.45077 - 0.05869 * (WL-WL_{max})]^{-1}]$$
This equation is derived from previous work showing a water level-depended
difference in sap flux between black ash and current co-dominant species
[@shannon-2018].

<!-- Run after set-up of inundation manuscript: 
sapmod <- 
  lm(saparea_cm2 ~ dbh_cm, data = read_csv("../Data/Sapwood_Area.csv"))

coef(sapmod)

inundation %>% 
  filter(site %in% c("140", "151", "152")) %>% 
  select(site, sampleDate, Species, pooledSpecies, sapFlux_g_cm2d1, rawWaterLevel_cm, transformedDz_kPa, dbh_cm) %>% 
  mutate(sapwood_area_cm2 = case_when(Species == "Yellow Birch" ~ 22.3 * dbh_cm - 236.1,
                                      Species == "Red Maple" ~ 17.04 * dbh_cm - 110.66,
                                      Species == "Black Ash" ~ 5.427 * dbh_cm - 43.65),
         sapflow_g_d1 = sapFlux_g_cm2d1 * sapwood_area_cm2) %>% 
  group_by(sampleDate, site, pooledSpecies) %>% 
  summarize(rawWaterLevel_cm = mean(rawWaterLevel_cm, na.rm = TRUE), 
            transformedDz_kPa = mean(transformedDz_kPa, na.rm = TRUE),
            sapflow_g_d1 = mean(sapFlux_g_cm2d1, na.rm = TRUE)) %>% 
  ungroup() %>% 
  mutate(normWaterLevel_cm = case_when(site == "140" ~ rawWaterLevel_cm - 11.2594,
                                       site == "151" ~ rawWaterLevel_cm - 1.51263,
                                       site == "152" ~ rawWaterLevel_cm - 6.81673)) %>% 
  pivot_wider(names_from = c(pooledSpecies), 
              values_from = sapflow_g_d1) %>%
  mutate(sapFlowRatio = `Non-Black Ash` / `Black Ash`) %>% 
  robustbase::glmrob(sapFlowRatio ~ normWaterLevel_cm,
                     data = .,
                     family = Gamma)
-->

Because individual models were trained for each wetland there are no training
data available for non-forested conditions under the field-study control sites.
Therefore we fit wetland water level models to only of 8 of the field-study
wetlands: the 8 treatment (girdle and ash cut) wetlands. Comparisons of the
combined and separate impacts of EAB and climate change were performed against
modeled baselines rather than field-study observed conditions. Limiting the
number of sites and comparing to modeled future simulated baselines with each
alternative vegetation condition ensures that comparisons are not biased by the
number of hydrology of sites within each group. Control conditions are
represented by the modeled baselines, avoiding the potential of identifying
model artifacts as significant when comparing observed and modeled results.

#### Future Climate Conditions

Wetland models were run under simulated current (1980-2009) and future
(2070-2099) conditions. Future conditions were evaluated under the RCP 4.5 and
RCP 8.5 representative concentration pathway scenarios [@vanvuuren-2011]. Daily
climate projections for all scenario-period combinations were taken from the
localized constructed analogs (LOCA) downscaling for the General Fluid Dynamics
Lab Coupled Model (GFDL-CM3) and National Center for Atmospheric Research
Community Earth System Model (CCSM4) global climate models (GCMS)
[@pierce-2014; @gent-2011; @griffies-2011]. LOCA downscaled data were retrieved
from the downscaled CMIP3 and CMIP5 Climate and Hydrology Projections archive
(http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/).

Performance of GCMs in the North American Great Lake region varies due to the
mostly unaccounted for regional climate impact of the Great Lakes
[@notaro-2015; @rood-2018]. GCM model selection was guided by performance and
sensitivity (magnitude of change under a given emissions scenario relative to
other GCMs) of the models in the Great Lakes region [@byun-2018]. Both GFDL-CM3
and CCSM4 were found to be well performing models for the Great Lakes region,
and CCSM4 and GFDL-CM3 represent a less and more sensitive models,
respectively. Inspection of LOCA daily simulations of GFDL-CM3 and CCSM4 for
the 1980-2009 period show excellent alignment with observed conditions during
the same period (Figure \@ref(fig:gcm-evaluation)). There appears to be fewer
dry days in the LOCA dataset relative to the observed dataset. Monthly and
seasonal precipitation totals showed good agreement between the LOCA and
observed datasets and no bias-correction was performed.

```{r gcm-evaluation, fig.cap = "Comparison of LOCA daily projections and observed values for the climate normal period (1980-2009) from the GFDL-CM3 (A-C) and CCSM4 (D-F) global climate models. Points show median value, bold and thin bar show 1 and 2 standard devations respectively, and filled polygons show distribution of values."}

knitr::include_graphics("output/figures/density_ridges_gcm_and_observed_climatology.tiff")

```

#### Stochastic Weather Generator  

The objective of projecting wetland conditions under future climate scenarios is
to understand the expected response, and the range of possible responses. The
30-year periods of daily projections provide a representation of the normal
climate variation under each scenario [@baddour-2007]. Quantifying the range of
weather patterns under future climate scenarios and the corresponding wetland
responses is best performed with many more years of simulation using stochastic
processes [@wilks-2012]. Stochastic weather generators (SWGs) are a tool for
generating synthetic time series of weather that simulate conditions under an
observed or projected climatology [@wilks-1999]. SWGs have also been used as a
downscaling technique for global or regional climate models [@verdin-2015;
@verdin-2019]. We followed the generalized linear model (GLM) based SWG
described in @verdin-2015, which is a type of the more general Richardson SWG 
[@richardson-1981].  

Our SWG consists of four GLM models to simulate precipitation occurrence,
precipitation amount, minimum daily temperature, and maximum daily temperature.
Minimum and maximum daily temperatures are fit using separate models with the
same model form. Each daily value is predicted using the previous day's minimum
and maximum temperatures, a pair of Fourier harmonics, and seasonal mean minimum
and maximum temperatures. To simulate the daily variation inherent within a
climate season additional error is added to the predicted daily minimum and
maximum temperatures. The additional error is drawn from a multivariate
skew-normal distribution conditioned on precipitation occurrence. This approach
captures the correlation between daily minimum and maximum temperatures and
between daily temperatures and precipitation occurrence. Precipitation occurrence
is modeled using the logistic family with a probit link function. Occurrence is
considered as a first order Markov process dependent upon the previous day's
precipitation occurrence [@wilks-1999]. Occurrence models also include a pair of
Fourier harmonics to capture any seasonality in precipitation patterns, and the
seasonal totals of precipitation for the observed period as covariates.
Precipitation amounts are drawn from a gamma distribution defined by shape and
scale parameters. The shape parameter is extracted from a gamma-family GLM with
a log link fit to the observed precipitation amounts. The predictors in the
model are a pair of Fourier harmonics and the seasonal precipitation totals. The
scale of gamma distribution describing the distribution of storm sizes is a
function of the GLM gamma-family shape and the predictors from GLM fit to 
precipitation amounts. This approach allows the
storm distribution to vary temporally and capture seasonality in not only total
precipitation but also the size of the storms. When seasonal means and totals
are included as predictors within the GLM predictors it is said that the SWG is
conditioned on the observed climate. When applying a conditioned SWG to a new
climate projection it is assumed that the distribution of the precipitation
amounts will remain constant under the new seasonally-defined climate scenario
[@wilks-2012]. We removed this assumption from our analysis by developing
separate SWGs for each model-scenario combination using the LOCA daily
projections for GLM fitting.

The work in @verdin-2015 shows that these SWGs can capture and simulate the
spatial correlation in temperature and precipitation. We did not include
spatial correlation in our SWGs because of the relatively small size of the
study area and the fact that each wetland is considered an independent system.
For each model-scenario combination the SWG was fit using LOCA daily
projections from GHCN-D sites USC00200718, USC00206220, USC00204104,
USC00476398, USC00207812, USC00201088, USC00204328, USC00475352, USC00206215,
USC00476398, USC00476518, USC00476939, USC00477140, USC00208680 to calculate
annual regional seasonal precipitation and temperature for conditioning data.
The projected daily data from the coordinates of the Bergland Dam
meteorological station were used to fit the SWG GLMs. Each SWG was used to
simulate 10,000 individual synthetic weather series with seasonal conditioning
drawn randomly from the 30 years of projection data. Conditioning the models on
individual years of observed data is intended to increase interannual (or
inter-simulation) variability, which can otherwise be limited in SWGs
[@wilks-1999].

<!-- 
It currently looks like the errors for the tmin and tmax models are too small 
when comparing the CDFs and PDFs of the SWG and LOCA tmin and tmax on a seasonal
basis. After some experimenting it looks like part of the problem might be that
the white noise added to the tmin/tmax simulations are from sd of monthly data,
which is smaller than sd of seasonal data. Using sd of seasonal data or
conditioning by month should solve this and get the KS tests to fail to reject
the NULL that the data are from the same distribution. Likely better errors 
could be computed using the skew_normal distribution. For JJA it looks like both
tmin and tmax have a negative skew that allows for almost perfect match with
PDF plot:
MASS::fitdistr(obs_data[gcm == "ccsm4" & scenario == "rcp45" & sample_season == "jja", tmin_c], dskew_normal, start = list(mu = 11, sigma = 4, alpha = -1), lower = list(mu = -Inf, sigma = .Machine$double.eps, alpha = -Inf))
The standard deviations of tmin and tmax may also vary with precip occurrence,
and SD of each of them would likely be related to each other. So the noise added
to the models should probably be mvrnorm() with a covariance matrix of tmin, tmax, 
and precip occurrence [@wilks-2012] 
-->

#### Data Analysis

Tests for statistical significance are unperformed and unreported with a single
exception. This research is entirely dependent on simulation data drawing
conclusions from modeled wetland water levels driven by synthetic weather
series generated from parametric descriptions of LOCA downscaled of GCM
projections. Reporting statistical significance and *p*-values would provide a
false sense of certainty and dichotomy to the future climate and wetland
conditions. In our results we report and account for the uncertainty inherent
in this type of simulation work. We draw conclusions by contrasting the
probability of observing certain wetland conditions under different vegetation
types and climate scenarios.

##### Total, EAB, and Climate Impacts

The distribution of modeled wetland water levels was evaluated for each future
vegetation and climate scenario combination. The 10,000 simulations for each
combination were summarized by simulated day of year to compute the median and
the 67% highest density intervals [@mcelreath-2020] for each combination. This
approach was used to quantify future conditions and baseline scenarios for
comparison and benchmarking. Modeled baselines provide consistency between the
baseline, or control, period and alternative vegetation cover and climate
scenarios, which is critical to drawing meaningful conclusions. <!-- TODO: find
citation about comparing future climate predictions to modeled baselines as
best practice --> Apart from the advantage of consistency, modeled baselines
can also provide a flexible tool to answer more questions about the drivers of
the the observed changes. Different baselines can be computed from the
simulations by combining the six vegetation-climate combinations. To evaluate
total impact of EAB and climate change, future non-black ash conditions were
compared to black ash under the current climate. Alternative baseline
comparisons include comparing each vegetation cover to itself under alternative
climate scenarios, and comparing vegetation covers to each other within a
climate scenario. These two alternative baselines allowed us to attribute how
much of the observed total impact was attributable to each driver.

##### Critical Ecohydrological Thresholds
Observed water levels within these wetlands show high interannual variation 
under field control and treatment conditions. Rather than consider this 
low-signal, high-variance variable, we defined critical ecohydrological 
thresholds (*CEHTs*) as a measure of impact. *CEHTs* were 
set to capture wetland **connectivity** to the downstream hydrologic network via 
streamflow and subsurface flow, **inundation** when wetland water levels were 
near the soil surface with surface water likely in microtopography hollows, and 
**drawdown** when wetland water levels dropped far below the surface of the 
wetland (Table \@ref(tab:ceht)). Wetland water levels 
were compared to these thresholds and the number of days that a wetland was 
above or a below a given threshold could be used to calculate the probability of
occurrence of that *CEHT*. We chose to calculate these probabilities at the 
monthly scale, though can than be computed for time scales from daily to 
annually.

```{r ceht, tab.id = "ceht"}

threshold_names <-
  factor(
    x = c("Connectivity", "Inundation", "Drawdown"),
    levels = c("Connectivity", "Inundation", "Drawdown"),
    ordered = TRUE
  )

water_levels <- 
  factor(
    x = c("≤ 5 cm below maximum site water level", "≤ 10 cm below wetland surface", "> 50 cm below wetland surface"),
    levels = c("≤ 5 cm below maximum site water level", "≤ 10 cm below wetland surface", "> 50 cm below wetland surface"),
    ordered = TRUE
  )

definitions <-
  factor(
    x = c("Wetland is contributing water downstream", "Wetland has surface flooding or soil saturation", "Wetland has dry surface and non-saturated soils"),
    levels = c("Wetland is contributing water downstream", "Wetland has surface flooding or soil saturation", "Wetland has dry surface and non-saturated soils"),
    ordered = TRUE
  )

data.frame(
  Name = threshold_names,
  Threshold = water_levels,
  Definition = definitions
) %>%
flextable::flextable() %>%
flextable::set_caption(
    caption = "Water level thresholds and defintions of three conditions considered as critical ecohydrological thresholds. The occurrence of these thresholds were used to compare the impacts of EAB and climate change on the future of black ash wetland hydrology.",
    autonum = officer::run_autonum(seq_id = "tab", bkm = "ceht")
  )

```


### Results

#### SWG Performance

Synthetic weather generator performance can be evaluated by comparing generated
weather sequences with the underlying observed/simulated climate data. The
summary statistics of the sequences should align well with underlying data and
the distribution of observed values should closely match the shape of the
underlying data [@gregory-1993]. Table \@ref(tab:swg-check-table) shows
seasonal and climate model scenario summaries of LOCA and SWG variables. Values
are reported as the mean and a 95% confidence interval (1.96 times standard
deviations) around the mean. The percent of outliers were calculated using the
intervals defined by the LOCA dataset. <!-- TODO: change mean ± sd to median
and HDCI --> When comparing all of the simulated SWG data (equivalent to 10,000
years) the mean values and interquartile ranges of daily minimum and maximum
temperature and precipitation align well with the LOCA-generated daily values
(Table \@ref(tab:swg-check-table)). Fewer total monthly precipitation outliers
were observed for the SWG synthetic series than the LOCA climatology. Minimum
and maximum temperatures showed an increase in outliers in the synthetic
weather series for JJA and SON seasons, and a decrease in outliers for the DJF
and MAM seasons. The difference in outliers was small and may be partially
explained by the scale of the two datasets. The LOCA data represents 30 years
of data compared to the 10,000 years of synthetic weather series. To compare
the SWG and LOCA values as climates, we randomly sampled 30 sets of 30
years from the observed 10,000 synthetic series. The distribution of these
samples were compared to the distribution of LOCA values in Figure
\@ref(fig:swg-check-plot). Figure \@ref(fig:swg-check-plot)A shows good
agreement for all seasons for daily minimum and maximum temperature. The
distribution of simulated monthly precipitation totals are shown in Figure
\@ref(fig:swg-check-plot)B. Relative to minimum and maximum temperature,
monthly precipitation totals show greater variability among the 30 sets of
30-year synthetic weather series. All three variables show that the synthetic
weather climates cluster around the LOCA value as a mean.

```{r swg-check-table}

targets::tar_read(swg_check_table) %>%
  flextable::set_caption(
    caption = "Comparison of 30 years of LOCA-generated climate values and 10,000 simulated weather series (SWG). Value are presented as mean with the 95% confidence interval in parentheses (± 1.96*SD). Percentage of outliers is the number of each value type that falls outside of the interval defined by ± 1.96*SD for the LOCA-generated values.",
    autonum = officer::run_autonum(seq_id = "tab", bkm = "swg-check-table")
  )

```


```{r swg-check-plot, fig.cap="Density plots of 30 years of LOCA generated data (bold), and 30 random samples of 30-year periods from the full 10,000 years of SWG synthetic weather series (light). Bold lines represent LOCA climatology, and fine lines show individual 30-year samples of SWG data."}

knitr::include_graphics("output/figures/density_lines_loca_and_swg_climatology.tiff")

```

#### Wetland Water Level Model Performance

Wetland water level model performance was evaluated using the retained
independent testing datasets. We calculated *R^2^* as a metric for the
relationship between observed and modeled values. *R^2^* cannot be used to
evaluate the accuracy of the models because consistent over or under predictions
can still result in high *R^2^* values [@krause-2005]. We use the median
prediction error to measure model bias and the root median squared error 
(RMedSE) to assess overall predictive accuracy. RMedSE was used in place of RMSE 
because we expect some outlier errors where regionally-informed rainfall records 
do not match observed wetland rainfall. In addition to the
overall accuracy of the models we want the models to accurately capture the
probability of occurrence of the *CEHT*s. We
tested the predicted probabilities of inundation, connectivity, and drawdown to
the observed probabilities of the same conditions.

One wetland, 119, had all of the withheld training years with *R^2^* below 0.6,
which was chosen as a threshold for unsatisfactory model results based on work
by @moriasi-2015 (Figure \@ref(fig:model-metrics)). Models at this wetland did
not show noticeably higher rates of bias or error relative to the other
wetlands. Site 119 is a closed basin with no surface outlet and was shown to
have different source water characteristics than other sites
[@vangrinsven-2017]. Table \@ref(tab:model-metrics-table) summarizes the wetland
model metrics by site status, presenting the minimum, maximum, mean, and median
values observed across all site-year combinations. Mean and median model
performance was the similar between site statuses for all model metrics
(\@ref(tab:model-metrics-table)). After excluding wetland 119 *R^2^* values ranged
from 0.44--0.95. The models had an overall negative bias, which indicates drier 
wetland conditions than were observed (Table \@ref(tab:model-metrics-table)). Both the 
median rRMedSE is moderate for both Control and Treated condition models (~10%),
but notable outliers were present (wetland 053) (Table \@ref(tab:model-metrics-table)).

```{r model-metrics, fig.cap="Wetland model performance metrics. Metrics were calculated for each site-year combination and are presented within sites. RMedSE and rRMedSE are root median squared error and root median squared error relative to wetland water level range within that year."}

knitr::include_graphics("output/figures/boxplot_wetland_model_metrics.tiff")

```
  

```{r model-metrics-table}

targets::tar_read(wetland_model_metrics_table) %>%
  flextable::autofit() %>%
  flextable::set_caption(
    caption = "Minimum, maximum, and median values of wetland water level model metrics, comparing observed and modeled wetland water levels. Metrics were evaluated on withheld test data for each combination of site and year of data. Prior to calculating these metrics three sites with incompatible hydrology for the model structure were removed from the analysis.",
    autonum = officer::run_autonum(seq_id = "tab", bkm = "model-metrics"))


```

When compared to observed wetland water levels, modeled wetland water levels
showed no significant difference in probability of occurrence of *CEHT*s (Table
\@ref(tab:predicted-probs-table)). Reported values and significant tests are from a linear
mixed model with dependent variable probability of occurrence, population
effects of observed/predicted and wetland status, and a group effect for
wetland [@lenth-2021]. Summary and statistical tests were performed using
estimated marginal means [@bates-2015]. The range of site-level probability of
each level of interest were similar between modeled and observed data (Figure
\@ref(fig:predicted-probs)). Though non-significant these results show
a slight systematic bias towards drier conditions in the modeled wetland water
levels. All remaining comparisons of current conditions to combinations of
future vegetation conditions and climate scenarios will use modeled current
conditions. This comparison ensures that the described non-significant model
bias does not impact results or conclusions. It also sets equal sample sizes
between the control and future conditions.

**Try vertical text for variable label then can horizontally expand LOCA and SWG values**

```{r predicted-probs-table}

targets::tar_read(wetland_model_predicted_probs_table) %>%
  flextable::autofit() %>%
  flextable::set_caption(
    caption = "Tests comparing predicted probability of occurrence for critical ecohydrological thresholds (*CEHT*). Connectivity, Inundation, and Drawdown respectively represent wetland water levels equal to or greater than modeled maximum water level, water levels within 10 cm of the wetland surface, and water levels more than 50 cm below the wetland surface. Reported values show the seasonal probability of daily water level exceeding the respective *CEHT*. Control and treated conditions refer to **observed** control and treated conditions or **predicted** modeled black ash and modeled non-forested conditions. All data presented derived from only the withheld test period for each wetland.",
    autonum = officer::run_autonum(seq_id = "tab", bkm = "predicted-probs"))

```

```{r predicted-probs, fig.cap="Tests comparing predicted probability of occurrence for critical ecohydrological thresholds (*CEHT*). Connectivity, Inundation, and Drawdown respectively represent wetland water levels equal to or greater than modeled maximum water level, water levels within 10 cm of the wetland surface, and water levels more than 50 cm below the wetland surface. Reported values show the seasonal probability of daily water level exceeding the respective *CEHT*. Control and treated conditions refer to **observed** control and treated conditions or **predicted** modeled black ash and modeled non-forested conditions. Black point and error bar represents estimated marginal mean and its 95% confidence interval. Individual points represents data observed or simulated for 1 year at 1 site. All data presented derived from only the withheld test period for each wetland."}

knitr::include_graphics("output/figures/jitterpoint_and_pointrange_wetland_model_predicted_probability.tiff")

```

#### Wetland Water Levels

Wetland water levels under future climate conditions were highly variable under
both scenarios and both vegetation types (Figure \@ref(fig:total-impact)).
Figure \@ref(fig:total-impact) shows the probability that wetland water levels
will be above or below the *CEHT*s defined above. The reported probability is
the proportion of observations where daily water level surpasses each threshold
out of all modeled daily water levels in a month (~240,000). The median
probability and a range representing 67% of all observations (highest density
continuous intervals, HDCI) are shown as point and line estimates for the less
and more sensitive future climate scenarios. Modeled black ash conditions are
also reported as the median and 67% HDCI, represented as a crossbar and shaded
area. Figures \@ref(fig:eab-impact) & \@ref(fig:climate-impact) present other
comparisons using the same symbology structure.

For non-forested conditions both future climate scenarios showed a decrease in
probability of connectivity relative to current wetland conditions (Figure
\@ref(fig:total-impact)A). The less sensitive climate scenario consistently
showed a lower median probability of connectivity and was more likely to have
simulated conditions outside of the HDCI of current conditions. The
more-sensitive climate scenario showed greatest agreement with current
conditions in July and August. Future forested conditions were more likely to
show connectivity in July and August under the more sensitive climate scenario,
and only slightly less likely than current conditions under the less sensitive
climate scenario (Figure \@ref(fig:total-impact)B). Connectivity to the larger
hydrologic network was much less likely to occur under the more sensitive
climate scenario and both vegetation conditions in May (Figure
\@ref(fig:total-impact)A & B). Under both vegetation conditions the
less-sensitive climate scenario resulted in a lower probability of connectivity
relative to the more-sensitive climate scenario. While both vegetation
conditions showed significant overlap between the less-sensitive climate
conditions and control conditions, the non-forested conditions tended to skew
towards current drier years (lower probability of connectivity) and future
forested conditions skewed towards current wetter years (high probability of
connectivity).

Inundation results in Figure \@ref(fig:total-impact)C & D closely resemble the
patterns observed in the connectivity results. Probability of inundation and
connectivity are not exclusive. This is inline with fit model parameters that
showed maximum sustained water levels were above the surface for almost every
wetland and treatment combination (**TABLE S1**). Generally, inundated
conditions would be less prevalent under future non-forested conditions than
the future-forested conditions. For non-forested conditions there was little
overlap between baseline conditions and the less sensitive climate scenario
results (Figure \@ref(fig:total-impact)D). This was contrary to non-forested
conditions where the less-sensitive climate scenario results showed more
overlap with baseline conditions, and the more sensitive scenario varied by
month relative to baseline conditions. Future-forested conditions showed
increased probability for inundation in July, August, and September and both
vegetation conditions showed lower probability for inundation in June (Figure
\@ref(fig:total-impact)D).

The probability that water levels will drop to less than 50 cm below the
wetland soil surface agree with the connectivity and inundation results (Figure
\@ref(fig:total-impact)E & F). Under both climate scenarios the future forest
conditions showed a lower probability of drawdown than under baseline
conditions, with almost no overlap with current conditions in July, August, and
September (Figure \@ref(fig:total-impact)F). Drawdown events were very unlikely
with near zero probability of occurring under future forested conditions under
both climate scenarios. Probability of drawdown in non-forested conditions
increased in magnitude and variability dramatically from July through August
under both future climate scenarios. Differences from current conditions were
greatest in August under the less sensitive climate scenario (Figure
\@ref(fig:total-impact)E).

```{r total-impact, fig.cap="Probability of observing water levels above (connectivity, inundation) or below (drawdown) ecohydrologically significant thresholds. Points represent median the median observed probability within a month and ranges represent 67% of the simulations as highest-density continuous interval (HDCI). Gray shaded area represents model simulations for control black (ash forest) conditions under the current climate and the black bar represents the control monthly median probability."}

knitr::include_graphics("output/figures/pointrange_and_slabinterval_ecohydrological_level_total_impact.tiff")

```

### Discussion

**Talk about more sensitive connectivity in May and inundation in June as indicators of earlier melt**
**Can tie it back to the sourcewater work for when snowmelt signal was present**

#### SWG

Our SWG performed well at simulating weather series that conformed to LOCA-downscaled future climate scenarios (Table \@ref(tab:swg-check-table), Figure 
\@ref(fig:swg-check-plot)). The synthetic weather series were well aligned with 
the LOCA series, with precipitation showing more inter-simulation variability 
than minimum or maximum temperature. This result is understandable when we 
consider that in the study region precipitation is only weakly seasonal, while 
minimum and maximum temperature have strong seasonal signals. <!--For the downscaled 
LOCA values the CV of the seasonal-mean precipitation was 0.24 cm and was 
3.45 and 0.75 &deg;C for the seasonal-mean minimum and maximum temperatures, 
respectively.--> There may be a slight trend in under estimating the range of 
observed future climate conditions (Table \@ref(tab:swg-check-table)). This is 
likely a result of reversion to the seasonal mean due to a lack of available 
external weather drivers. Under future climate scenarios extreme events are 
expected to be more commonplace (**CITATION**), and these results suggest event size 
could be better parameterized by increasing the probability of extreme events 
through selection of distribution choices and modeling framework [@verdin-2019]. 
This shortcoming of the SWGs should not affect our findings because our
objective is not to model response to extreme events, but overall probability of
*CEHT* water levels. These water levels will be affected by extreme
precipitation events but the systems do not have infinite storage and more
extreme events can be expected to pulse through the system more quickly
(**CITATION**). The quicker the pulse passes through the system the smaller the
impact on overall probability of wetland water levels exceeding *CEHTs*

<!-- Get the CVs
loca_simulations[
  i = station_name == "bergland_dam",
  by = .(gcm, scenario, climate_season = as.climate_season(sample_date, TRUE)),
  j = lapply(.SD, function(x) {mean(x)}),
  .SDcols = c("precip_cm", "tmax_c", "tmin_c")
  ][
  (gcm == "ccsm4" & scenario == "rcp45") |
    (gcm == "gfdl-cm3" & scenario == "rcp85"),
  j = lapply(.SD, function(x) {sd(x) / mean(x)}),
  .SDcols = c("precip_cm", "tmax_c", "tmin_c")
]
-->

We chose to use an SWG approach that relied on only a single estimate of climate 
within the study area. Alternatively, regional variation in climate and 
weather patterns could be incorporated by using additional meteorological 
stations in the development of the SWG. A spatially-explicit SWG approach 
takes into account the covariance of regional meteorological stations and has 
been shown to faithfully capture larger regional trends in climate and weather 
[@verdin-2015]. Within our study region there are no clear weather or climate 
gradients and so we assumed that each study wetland was within a single 
climatic unit. This results in a much simpler SWG-construction process, but the 
resulting SWG cannot be used outside of the study region first refitting on data
from a different meteorological station(s). To simulate black ash wetland 
response on a regional basis would require addressing known temperature and 
precipitation gradients through the range of black ash. 

In addition to using only a single spatial representation of future climate , we
chose to use individual years of simulation. The direct impact of this is to
missing long-term droughts and wet spells. These events are expected to increase
in frequency under future climate scenarios (**CITATION**). We deliberately
chose to use single-year simulations because of the length of available training
data. During our study the region was in the process of transitioning from a
relatively dry period to a relatively wetter period (\@ref{fig:water-balacne}).
Our training data are primarily drawn from the drier portion of the study
period. Without longer-periods of training data we do not feel that we can
capture both the intra- and inter-annual dynamics of wetland hydrology. If we
assumed all 10000 years of synthetic weather series were contiguous weather
records, the longest period of drought years (annual precipitation less than
annual PET) was 7 and 9 years for the less and more sensitive climate scenarios,
respectively. By using single-year simulations we are implicitly assuming that
fall and winter (SON and DJF) precipitation is high enough to recharge the local
hydrology driving these wetlands, restarting the cycle at or near maximum
wetland water levels. When we consider only fall and winter precipitation we
found 0%, 8.7%, and 10.1% of synthetic weather series had totals below the tenth
percentile of observed data, less-sensitive, and more-sensitive climate
scenarios, respectively. Unsurprisingly close to 10% of the synthetic weather
winters were drier than the tenth percentile of the LOCA downscaled values.
Synthetic (and LOCA) winters were wetter than observed winters, supporting the
assumption that dormant-season precipitation would be enough to recharge wetland
water levels.

<!-- Drought Years
swg_simulations_loca[
  (gcm == "ccsm4" & scenario == "rcp45") |
    (gcm == "gfdl-cm3" & scenario == "rcp85"),
  .(drought = sum(precip_cm) < sum(abs(pet_cm))),
  by = .(gcm, scenario, simulation_id)
][,
  .(max_drought = {
    rles <- rle(.SD[["drought"]])
    max(rles$lengths[rles$values])
    }),
  by = .(gcm, scenario)]

-->

<!-- Percentage of dry winters simulated compared to loca simulation and observed

obs_dry <- external_met[, .(precip_cm = sum(precip_cm)),
             by = .(sample_year, station_name, climate_season = as.climate_season(sample_date, TRUE))
             ][, .(obs_dry_precip_cm = quantile(precip_cm, probs = 0.1)),
               by = .(climate_season)]

loca_dry <- 
  loca_simulations[
    ((gcm == "ccsm4" & scenario == "rcp45") |
      (gcm == "gfdl-cm3" & scenario == "rcp85")) & 
      station_name == "bergland_dam",
    .(precip_cm = sum(precip_cm)),
    by = .(gcm, scenario, sample_year, climate_season = as.climate_season(sample_month, TRUE))
  ][, .(obs_dry_precip_cm = quantile(precip_cm, probs = 0.1)),
    by = .(gcm, scenario, climate_season)]

swg_simulations_loca[
  (gcm == "ccsm4" & scenario == "rcp45") |
    (gcm == "gfdl-cm3" & scenario == "rcp85"),
  .(precip_cm = sum(precip_cm)),
  by = .(gcm, scenario, simulation_id, climate_season = simulation_season)
  ][
    loca_dry,
    on = c("gcm", "scenario", "climate_season")
  ][climate_season %in% c("son", "djf"),
    .(precip_cm = mean(obs_dry_precip_cm > precip_cm)),
    by = .(gcm, scenario)]

swg_simulations_loca[
  (gcm == "ccsm4" & scenario == "rcp45") |
    (gcm == "gfdl-cm3" & scenario == "rcp85"),
  .(precip_cm = sum(precip_cm)),
  by = .(gcm, scenario, simulation_id, climate_season = simulation_season)
][
  obs_dry,
  on = c("climate_season")
][climate_season %in% c("son", "djf"),
  .(precip_cm = mean(obs_dry_precip_cm > precip_cm))]
-->

```{r water-balance, fig.width=6, fig.height = 3, fig.cap = "Regional cumulative water availability for the period from November 1, 2005 through October 21, 2020. Water availability was calculated as the difference between rainfall plus melt and potential evapotranspiration."}

knitr::include_graphics("output/figures/line_plot_study_period_water_availability.tiff")

```

#### Ecosystem Specific Yield

The approach used to calculate *E~Sy~*, relating the magnitude of
water level change to the drivers of that change was an alternative to those
presented in previous research. We developed this alternative based on
unsuccessful attempts to implement *E~Sy~* equations through the rainfall/rise
method described in McLaughlin, et al (-@mclaughlin-2014). Directly proving the
theoretical basis of this approach is outside the scope of this study, and would be best
performed with specifically-designed experimental or simulation work. Our research has
however provided empirical evidence for the application of our approach. Section
CASE_STUDY in Chapter TRANSDUCER_ERRORS demonstrates that *ET* estimates derived
using this approach perform comparably to other implementations for calculating
*E~Sy~* (see White method details above). Those results showed that calculated
*ET* was correlated with daily *PET* estimates from nearby meteorological
stations. To extend that comparison we examine how well a simple regression between daily
water level change and adjusted or unadjusted hydrologic drivers (*PET*, *R*, *M*) perform. Figure
\@ref(fig:esy-impact) shows the results of a quantile regression of pooled data,
(not used elsewhere in this study) explaining daily change in wetland water level as a
function of *PET* *P* and *M*. We found that *pseudo-R<sup>2</sup>* of the model
improved from 0.53 to 0.63. The difference between *pseudo-R<sup>2</sup>* 
are large enough to be of note, but more importantly we can see that water level drawdown is
extremely underestimated by the unadjusted model (Figure
\@ref(fig:esy-impact)A). Our results show that periods with the highest drawdown
are not correctly modeled using the raw inputs. These periods of high drawdown
are predicted when the *E~Sy~* adjusted hydrologic inputs are used. Periods of
high drawdown generally occurred mid-season when water levels were farther below
the surface. Our approach agrees with others [@mclaughlin-2014; @watras-2017] finding that the
impact of *E~Sy~* increases as water levels drawdown and that flooded conditions
approach an *E~Sy~* of 1 [@mclaughlin-2014]. This relationship results in the largest
adjustment to inputs during these mid-season periods, explaining the differences
seen between panels A and B in Figure \@ref(fig:esy-impact). These results
together with those in Section CASE_STUDY in Chapter TRANSDUCER_ERRORS provide
empirical evidence that this alternative formulation captures the variation of
*E~Sy~* with wetlands and can be used in scenarios where the rainfall-runoff
ratio approach is masked by outside influences such as low-frequency hydrologic
signals (@zhu-2011) and surface water connectivity (@watras-2017).

<!-- TODO: What if the response of black ash transpiration to lowering water levels (rapid 
increase in T) is actually the signal I am capturing, not *E~Sy~*? This may be 
partly happening, but the response between rainfall and water level response is 
well explained by the current *E~Sy~* response, suggesting that the modeled 
phenomenon is primarily *E~Sy~*.-->

```{r esy-impact, fig.width=6, fig.height=8, fig.cap="Predictions of daily change in water level using a a quantile regression with observed daily water level as the dependent variable and *PET*, *rainfall*, and *snowmelt* as the independent variables. The model was fit two times, once with raw hydrologic inputs (A) and once with inputs adjusted by ecosystem specific yield (B). The right panels show the full dataset while the left panels zoom in on a interval [-5, 5] for both the observed and predicted water levels. The dashed line shows a 1:1 line."}

knitr::include_graphics("output/figures/scatterplot_and_zoom_delta_water_level_esy_predictions.tiff")

```


#### Wetland Water Level Models

<!-- 
wetland_model_metrics[
  metric == "r2" & !(site %in% c("119")),
  mean(value > 0.6)]
-->

Wetland model performance was adequate for the objective of this study, which
was to observe changes in the hydrographs of black ash wetlands and compare the
probability of occurrence for *CEHT*s. Overall the retained models showed good
correlation between observed and modeled water levels (Figure
\@ref(fig:model-metrics)A, Table \@ref(tab:model-metrics-table)). We found that 91% of
site-year combinations outside of site 119 had *R^2* values above 0.6 for the
model test period. RMedSE and relative RMedSE had similar patterns among sites (Figure
\@ref(fig:model-metrics)). This suggests that from both absolute and relative
standpoint, the wetland models as developed capture dynamic wetland water levels
more accurately than static wetland water levels. Daily percent bias estimates
(PBIAS) would be helpful to evaluate error relative to observed daily water
level. However, wetland water levels cross 0 so that the relative error
approaches infinity as water levels rise or fall towards the surface.

Importantly for our conclusions the wetland models performed well predicting the
probability of occurrence of *CEHT* levels. Table \@ref(tab:predicted-probs-table) and
Figure \@ref(fig:predicted-probs) show how probability of occurrence of
connectivity, inundation, and drawdown were not found to be different between
the observed and modeled data in the test periods. They also show that there is
no non-significant systematic model bias (positive or negative) in the
probabilities. Although no systematic model bias was identified for *CEHT*
analysis, Table \@ref(tab:model-metrics-table) does indicate a bias towards drier
simulations when daily water level measurements are compared. Therefore the use
a modeled control baseline is still an important safeguard against drawing
conclusions based on potential from model artefacts. Comparing the
systematically-biased results could exaggerate or mask real expectation of
wetter or drier wetland conditions.

The next step for these models is to fit them using a hierarchical Bayesian
approach. This should improve model fit by incorporating information from all
sites into a single model for each vegetation condition while allowing model
parameters to vary between sites. This approach will also increase the available
data for each vegetation condition when fitting the models as well as the
overall available data. With three hierarchical models, the field-study control
sites can incorporated without resulting in unbalanced vegetation-conditions
comparisons. If model performance does not significantly improve with that
approach more, or potentially all years of data can be used to parameterize the
wetland models, rather than the current single year of data for each period.
Using most or all of the data for model parameterization precludes the options
of evaluating model performance by cross-validation. Kozak and Kozak
(-@kozak-2003) provide an alternative view that replaces cross-validation with
lack-of-fit metrics, which evaluate model predictions for violation of model
assumptions (e.g., non-random error). Their main conclusion was that
cross-validation increased the variance of model estimates potentially resulting
in poor model choice.

#### Future Hydrologic Conditions

The importance of these simulation results is to help relate potential future
simulations to current conditions in dry, wet, or 'normal' years. This allows
researchers and managers to answer questions similar to: What will these sites
look like in a 'normal' year in the future? How will these sites respond to
wet/dry conditions in the future? For example, we could expect that during a wet
year under the less sensitive future climate scenario a non-forested wetland on
these sites would have a similar probability of surface inundation as a black
ash wetland does today (Figure \@ref(fig:total-impact)C). In general, black ash
wetlands that become non-forested can be expected to be have wetter conditions
under more sensitive future climate scenarios, and under less-sensitive future
climate scenarios are more likely to experience conditions similar to today.
Conclusions could be reached for each combination of climate and vegetation
scenario and this information can and should be used to influence management
approaches and timing in responding to EAB in black ash wetlands.

The combined impact of EAB and climate is shown as the difference between future
climate conditions (blue and orange) and black ash conditions under historical
climate (shaded gray) in Figure \@ref(fig:total-impact). We can draw additional
information about the individual effects of future climate and vegetative cover
on wetlands by altering how we calculate the baseline. In Figure
\@ref(fig:total-impact) the baseline is set as simulated black ash forests under
current climate conditions. To isolate the effects of each climate scenario of
wetland hydrology we can compare the response of a vegetative under future
climate conditions against its response under current climate conditions.
Conversely we can evaluate the impact of vegetative cover by comparing each
vegetative cover to simulated black ash forests under current climate
conditions. Figures \@ref(fig:climate-impact) and \@ref(fig:eab-impact) show the
probability of occurrence of each *CEHT* under conditions where the vegetative
cover and climate, respectively, are held constant. Table
\@ref(tab:impact-attribution) presents the change in probability of occurrence
for each critical water level and alternative vegetative condition that can be
attributed to EAB or a less or more sensitive future climate scenario.

We found that EAB impacts lead to slightly wetter conditions under non-forested
conditions and the present climate. This agrees with observed water level
response analysis in Van Grinsven, et al (-@vangrinsven-2017), where absolute
water level response to simulated post-EAB conditions resulted in wetter
conditions and significantly lower growing season drawdown rates. A modeled
alternate forest composition (Future Forested) shows that EAB impact alone would
lead to dramatically wetter conditions relative to black ash forests. This
result is expected as we modeled our future forested conditions based on
previous work showing that co-dominant hardwoods had significantly lower
seasonal transpiration estimates than black ash (@shannon-2017). Under both
alternative vegetative conditions we found that the isolated impact of future
climate scenarios would lead to much drier sites. The probability of surface
water connectivity or inundation occurring are dramatically decreased under both
climate scenarios. The probability of wetland water levels dropping below 50 cm
were increased under future climate scenarios for the non-forested conditions,
but showed little change for non-forested conditions. The probability of each 
*CEHT* is not symmetrical around the median, indicating for instance that the 
impact of EAB on Connectivity is more likely to be wetter conditions than 
looking at the median alone would suggest (\@ref(tab:impact-attribution), 
\@ref(fig:eab-impact)).

```{r climate-impact, fig.width=6, fig.height=3, fig.cap="Probability of occurrence for critical ecohydrological wetland water levels under combinations of vegetative cover and climate conditions. Each column compares simulated wetland response within a fixed vegetation condition across multiple climate scenarios. Probability of occurence is the proportion of days in each month simulated water levels (10000 simulations) reached each threshold water level. Connectivity: connected to surface drainage; Inundation: within 10 cm of soil surface; Drawdown: more than 50 cm below soil surface. Values are reported as the median and the bounds of a 67% highest density continuous interval (an interval, potentially asymmetric, that contains 67% of simulated values)."}

knitr::include_graphics("output/figures/pointrange_and_slabinterval_ecohydrological_level_climate_impact.tiff")

```

```{r eab-impact, fig.width=6, fig.height=3, fig.cap="Probability of occurrence for critical ecohydrological wetland water levels under combinations of vegetative cover and climate conditions. Each column compares simulated wetland response within a fixed climate scenario across multiple potential vegetative covers. Differences within a column represent the impact of EAB and potential management decisions. Probability of occurence is the proportion of days in each month simulated water levels (10000 simulations) reached each threshold water level. Connectivity: connected to surface drainage; Inundation: within 10 cm of soil surface; Drawdown: more than 50 cm below soil surface. Values are reported as the median and the bounds of a 67% highest density continuous interval (an interval, potentially asymmetric, that contains 67% of simulated values)."}

knitr::include_graphics("output/figures/pointrange_and_slabinterval_ecohydrological_level_eab_impact.tiff")

```

```{r impact-attribution}

targets::tar_read(impact_attribution_table) %>%
  flextable::set_caption(
    caption = "Expected change in probability of occurrence for critical ecohydrological water levels that can be attributed to future climate scenarios or emerald ash borer (EAB) impact and management response to EAB. EAB impact values represent the difference between simulated black ash forests and each alternative vegetative condition under current climate conditions. Climate impacts respresent the difference between each alternative vegetative condition under the current climate and future climate (less or more senstive) scenarios. Differences are calculated using the daily values from July through September, highlighting growing season impacts. Connectivity: connected to surface drainage; Inundation: within 10 cm of soil surface; Drawdown: more than 50 cm below soil surface. Values are reported as the median and the bounds of a 67% highest density continuous interval (an interval, potentially asymmetric, that contains 67% of simulated values).",
    autonum = officer::run_autonum(seq_id = "tab", bkm = "impact-attribution"))

```

Our hypothesis that the impacts of EAB and climate change would counteract each
other is partially supported by our results. The dynamics of wetland water
levels are complex and climate change and invasive species can each have major
consequences [@burkett-2000; @moomaw-2018]. We can summarize our findings by
saying we expect wetlands that transition from black ash forests to non-forested
sites to become drier, primarily due to the effects of increased summer
evaporative demand (Figures \@ref(fig:total-impact) & \@ref(fig:future-climate),
Table \@ref(tab:impact-attribution)). Our simulations saw less summer
precipitation under the less sensitive scenario, while under the more sensitive
scenario we observed an increase in large storms (Figure
\@ref(fig:future-climate). If management or natural regeneration leads to
establishment of future forests with lower site transpiration rates, these
wetlands can be expected to have the occurrence of high water conditions remain
stable or increase in frequency, while dry conditions (water levels more than 50
cm below the surface) become more rare (Figures \@ref(fig:total-impact) &
\@ref(fig:future-climate), Table \@ref(tab:impact-attribution)). Our results
show that the impact of climate change in this region will lead to consistently
drier wetland conditions. Changes to vegetative cover can amplify or counteract
hydrologic alterations attributable to climate change.

```{r future-climate, fig.width=6, fig.height=3, fig.cap="Distribution of 30-day total precipitation and potential evapotranspiration presented by climate season and climate scenario. All data present is fro the Bergland Dam, MI, USA GHCND station. Observed climate is from 1980-2009, future climate is from LOCA-downscaled climate scenarios for the period from 2070-2099. Potential evapostranspiration was calculated from observed and GCM outputs using the Hargreaves-Samani equation ([@hargreaves-2003])."}

knitr::include_graphics("output/figures/density_ridges_future_and_observed_climatology.tiff")

```


<!--
The high interannual variability of these systems does not show much
differentiation between climate scenarios and vegetation conditions. Examining
the probability of occurrence of certain threshold water levels does show
differentiation between climate scenarios and vegetation conditions. Similar to 
how interannual differences and hydrologic drivers masked treatment effects in 
raw water level signal.

Wetland water levels under observed conditions were highly variable
[@vangrinsven-2017] similar to simulations.

Discuss skew of distributions. Some months (check panel C in probs figure) exhibit more positive or negative 
skew in the observed probabilities.

Observed water levels are highly variable so it's not surprising that future 
conditions are also highly variable.

Connectivity of landscape wetlands [@creed-2014]
-->

#### Drivers of Future Wetland Hydrologic Conditions

##### Future Climate

The future climate scenarios used in this study present two alternative future
scenarios. The less sensitive scenario summer is projected to be slightly drier
than observed conditions, with similar *PET* (Table \@ref(tab:future-climate-table).
However, when summer water deficit is calculated as the *P-PET* we see that with
a given year conditions will show a decrease in water availability. The more
sensitive scenario is a much wetter projection than both the observed climate
and future climate (Table \@ref(tab:future-climate-table). There is little to no
overlap in the 67% HDCIs for summer precipitation between the more sensitive and
other scenarios. A commensurate increase in summer evaporative demand may lead
to slightly wetter summer conditions under the more sensitive scenario. Spring 
and early summer water levels are currently strongly influenced by the snowfall 
and melt regimes [@vangrinsven-2017]. Our simulations show that snowmelt was 
reduced in quantity and shifted earlier in the year (Supplemental Figure S1) for
both future scenarios, with more dramatic impacts under the more sensitive 
scenario. This shift will likely have the result of increasing water deficit on
the landscape relative to the values presented in Table \@ref(tab:future-climate-table).

```{r future-climate-table}

targets::tar_read(future_climate_table) %>%
flextable::autofit() %>% 
flextable::set_caption(
    caption = "Summer (JJA) precipitation, potential evapotranspiration, and water deficit (P - PET) for the climate scenarios in this study. Future climate data were taken from LOCA downscaled daily values, and observed climate was taken from observed conditions at Bergland Dam, MI, USA from 1980-2009. Reported values are the median and bounds of the 67% highest density continous interval of seasonal total precipitation by observed or projection year.",
    autonum = officer::run_autonum(seq_id = "tab", bkm = "future-climate")
  )


```

##### Interaction of *AET* and *E~Sy~*
It is expected that vegetation conditions with lower transpiration rates have
lower rates of water level drawdown and therefore lower increased water levels.
In these wetlands this effect is compounded by the dynamics of *E~Sy~*. *E~Sy~*
increases in magnitude as water levels decline (Figure \@ref(fig:esy-function)),
which may be a result of reduced soil pore space and wetland cross-sectional
area so that a smaller volume of change results in a large water level
difference [@mclaughlin-2014]. The effect is that *PET* results in smaller changes
to wetland water levels under wet conditions than dry conditions. In this way
the impact of EAB may result in to a feedback loop in which water levels remain
elevated:

1. *AET* begins to draw down water levels, *E~Sy~* is low,  
2. Water levels decline slowly because of reduced *AET* and low *E~Sy~*,  
3. As *PET* increases to mid-season peak, *E~Sy~* and *AET* are both lower than under black ash,  
4. High *PET* impact on water levels are reduced by lower *AET* and sustained 
low *E~Sy~*  

Under the current black ash conditions *E~Sy~* increases more rapidly due to
higher black ash *AET*. The increase in *E~Sy~* accelerates the impact of higher
*AET*, creating faster and larger declines in wetland water level throughout the
growing season (Figure \@ref(fig:pet-impact)). This feedback loop may partially
explain the persistence of hydrologic impact following EAB disturbance observed
in Michigan and Minnesota [@diamond-2018; @kolka-2018]. In effect, the impact of
EAB may shock these systems into an alternative stable state of elevated water
levels [@scheffer-2003]. Our results indicate that future climate scenarios will
likely have a large enough impact to again shock these systems out of a stable
state. We cannot say from these simulations whether the wetlands will reach a
new stable state under future climate conditions or what that state would be.

```{r pet-impact, fig.cap="Modeled and effective actual evapotranspiration under various climate scenarios and vegetation conditions. Modeled AET is the synthetic PET multiplied by the model coefficient used in the relevant wetland model. Effective AET is the modeled AET times the value of *E~Sy~* predicted from the contemporaneous water level. The difference between the two values demonstrates the two-stage impact of reduced ET, where AET and E~Sy~ are reduced."}

knitr::include_graphics("output/figures/lines_pet_demand_by_cover_and_climate.tiff")

```

##### Reduced Evaporation and Non-Canopy Transpiration

We observed wetter conditions observed under our Future Forest simulations
relative to both Black Ash and Non-Forested simulations under all climate
scenarios (Figures \@ref(fig:predicted-probs), \@ref(fig:eab-impact)).
Two factors are likely contributing to this result. The first is that our future
forest composition is known to have lower transpiration rates than the existing
black ash canopy [@shannon-2018]. As a model choice we opted to assume that the
next most dominant canopy species would be the likely replacement canopy
species. The future forest we modeled is a narrow range of potential future
forest compositions. natural regeneration or planting efforts may lead to future
forested species composition that more closely matches the evaporative potential
of the current black ash canopy. Secondly, the future forested conditions may
have a lower total *AET* than the non-forested conditions. In black ash wetlands
in Minnesota simulated post-EAB conditions (girdled and standing ash) were found
to have higher water levels than site were the ash stems were harvested and
removed [@diamond-2018]. The authors attributed the result to limited *AET*
leading to reduced solar energy and wind-driven boundary layer mixing due to the
still standing stems. Our conditions have a notable difference from that study
in that we are assuming there are living trees on site. However a closed canopy
would have the same effect of reducing open water evaporation and understory
transpiration. With sufficiently reduced canopy transpiration the effect of
certain forest compositions could result in reduced overall *AET*. This was
observed in reverse in northern Minnesota when total *AET* increased after
low-transpiration black spruce was removed from a peatland (**CITATION**).

The contrast between the Non-Forested and Future-Forested simulated hydrology
suggests an important management tool. Iverson et al [-@iverson-2016] and the
work from Bolton [-@bolton-2018] and Looney [-@looney-2015] laid out a framework
and results for evaluating potential replacement species considering site
conditions of black ash wetlands. The results presented here show the
opportunity for management decisions that take into account the impact of future
vegetation on site hydrologic conditions. The general trend of drier future
conditions can be to some degree counteracted by management for cover with lower
evapotranspiration rates. Drier conditions are not always preferable across the
landscape and this tactic could be used to retain water on the landscape,
creating refugia of standing water or cool moist soils for flora and fauna.
Although much of the region is expected to have drier summers our localized
study area is expected to have wetter summers. Our results can still be applied
to the broader region. The impact of climate change in other areas of the region
would become even more pronounced. Vegetation conditions would still amplify or
counteract climate impacts, and in some areas may become even more important to 
providing wet or moist refugia on the landscape.

### Future Research

It should be stated that these models and simulations cannot capture the full
range of potential changes from EAB, other potential invasive species, and
climate change. The following is a non-exhaustive list of processes and assumptions 
that this work could not address and warrant further study.

- Watras [-@watras-2014] found that there are decadal scale oscillations of water levels 
within the Great Lakes and inland lakes in the Great Lakes Region. These can be 
expected to impact wetland water levels similarly.  
- While our results showed a smaller and earlier snowmelt pulses
we did not attempt to model the more complex dynamics of rain-on-snow events and 
a compressed melting period reducing recharge to the intermediate groundwater 
sources shown to feed these wetlands [@vangrinsven-2017].  <!--TODO: snowmelt and ground citations--> 
- Our single-year simulations do not account for the potential effects of 
multi-year drought or for periods of significantly reduced fall/winter 
precipitation on wetland water level rebound.  
- The response of each species and plant community will show individual climate
change responses that are not captured in these models and may be non-linear
[@short-2016].  
- Beyond the potential for unforeseen vegetation responses to future climates is
the uncertainty of future climate conditions. The authors of the fourth national
climate assessment have highlighted that early climate prediction have
under-predicted contemporary shifts in response to climate change (TODO: add
citation from appropriate chapter). Unfortunately this means that our models
built on those simulations may be underestimating the magnitude of future
changes.  

Wetlands provide a range of ecological services from water quality and stormflow
retention to carbon storage, all of which stand to be impacted by climate change
[@moomaw-2018]. We did not attempt to quantify how changes in wetland hydrology 
will influence these other services. By highlighting the magnitude and drivers 
of hydrologic change in these systems we hope to spur future research.

### Conclusions

<!-- TODO: reference these from intro "3) site conditions following EAB infestation, and 4) site conditions in a future climate." -->

Our research has shown that changes in evaporative demand and precipitation
regimes will likely result in drier conditions on what are now black ash
wetlands. In addition the functional loss of black ash due to EAB presents 
challenges and opportunities for the future of these wetlands. A large extent of 
forested wetlands may require intervention to retain desired benefits or 
features they currently provide. These interventions can be used to drive sites 
towards wetter conditions, retaining water on the landscape that may otherwise 
be lost under a changing climate.  

To implement long-term management objectives in wetland and forest management 
requires a deep understanding of the systems. We have shown that as water levels 
decline, the impact of each additional driver increases due to changes in
*E~Sy~*. Changes in *AET* can interact with *E~Sy~* to create feedback loops 
underscoring that knowledge of species adaptation to wet conditions and capacity 
to respond quickly to drying conditions is critical for projected future site 
conditions. The transition from black ash to alternative vegetative cover can 
counteract or amplify the impacts of future conditions with higher evaporative 
demand. And the relationship between evaporative demand and *E~Sy~* can amplify 
the effect of wetland water level drawdown.

### Acknowledgements
I would like to thank Andrew Verdin for discussion of synthetic weather 
generators and code examples. Additionally my wife Danielle Shannon for both the 
support she provided during writing, and for share her knowledge of evaluating 
and communicating climate change impacts on forests.

### References