Forecasting the flooding dynamics of flatwoods salamander breeding wetlands under future climate change scenarios

Study sites

We conducted our work on Eglin Air Force Base (Eglin, FL, USA) in the Florida panhandle. Access to field sites was approved by the US Fish and Wildlife Service and Jackson Guard (Eglin’s Natural Resources Division; Cooperative Agreement Number F14AC00068). Eglin is a large military installation (>187,000 ha) located in Okaloosa, Walton, and Santa Rosa Counties and consists primarily of actively managed longleaf pine forests. The wetlands studied here have been the subject of ongoing research to understand the variation in hydrologic regime (e.g., Chandler et al., 2016; Chandler et al., 2017) across the landscape with the goal of managing and conserving valuable breeding habitat for flatwoods salamanders (see Chandler et al., 2017 for additional study site details). Extensive habitat management in both uplands and wetlands has made Eglin one of the few remaining strongholds for flatwoods salamanders (Gorman, Haas & Himes, 2013; U.S. Fish and Wildlife Service, 2020).

Water level data

As part of our long-term work on flatwoods salamanders, we previously installed water level monitoring wells at the approximate deepest point in 35 wetland basins on Eglin. All wetlands were located within 20 km of one another, and wetland areas ranged from 0.19–20.92 ha, although all but two wetlands were less than six hectares. Wells were installed from November 2014 through December 2017 such that we had between 1,061 and 2,338 daily measures of water level for each well (see Table S1 for more details). Each well consisted of a 3.8 diameter screened PVC pipe 1 m below ground with a HOBO U20 pressure transducer (Onset Computer Corporation, Bourne, MA) at the bottom of each well that recorded total pressure and temperature at 15-minute intervals. Total pressure data were corrected for barometric pressure variation using pressure sensors installed following methods in McLaughlin & Cohen (2011) at two locations within 9 km of any given wetland (Table S1). For each wetland, we used the resulting 15-min water level data to calculate mean daily water level (mm; positive indicating water levels above ground surface), along with the timing and duration of flooded conditions.

Current climate data

To determine precipitation inputs, we installed four rain gauges (HOBO Data Logging Rain Gauge RG3-M) such that the distance between any of the 35 wells and a rain gauge was no more than 4 km (Table S1). Three of the rain gauges were installed on 20 November 2014 and the fourth was installed on 19 June 2018. We used data from these gauges to calculate daily rainfall (mm), cumulative rainfall over the previous seven days (mm), and the number of days since the last rain event.

We obtained daily minimum and maximum temperature (°C) data from the Florida Automated Weather Network (FAWN) over the course of our study period using the Jay weather station (Station ID 110), which is located 32–58 km from each wetland. We calculated daily potential evapotranspiration (PET, mm) from minimum and maximum temperatures using the Hargreaves-Samani method (Hargreaves & Samani, 1985). Lastly, we calculated the 12-month standardized precipitation-evapotranspiration index (SPEI) for each well from 1 January 1981 to 31 May 2022. The SPEI is useful for assessing drought severity (WMO, 2006) and quantifying and comparing water balances across locations (Stagge et al., 2014). To calculate monthly SPEI, we used daily minimum and maximum temperature data obtained from PRISM Climate group (Oregon State University, http://prism.oregonstate.edu, created 26 Oct 2022), which were necessary because the FAWN data do not have a long enough time series for the SPEI calculation. For the same time period, we used the Hargreaves method (Hargreaves, 1994) to calculate monthly PET, which was used for calculating SPEI.

Hydrologic model

We developed a statistical model to predict daily water levels at each of the 35 wells using a Bayesian first order autoregressive fixed effect model, which was based on visual inspection of partial autocorrelation plots for each well. For each wetland basin, the water level of a given day was modeled as a function of the first order autoregressive term (AR), the amount of rain received during the previous day (Precip), and their interaction (Precip:AR). Models also included linear and quadratic terms for PET, the sum of the rain received in the previous seven days (WeekPrecip), and SPEI. We also included interactions between SPEI and the amount of rain received during the previous day (SPEI:Precip), and the interaction between SPEI, the first order autoregressive term, and the amount of rain received during the previous day (SPEI:AR:Precip). Therefore, for each of the 35 well-specific models, we estimated 11 parameters, including an error term, as shown in the following equation. (1)${\displaystyle {\text{WaterLevel}}_t=\alpha +\beta AR+\beta Precip+\beta PET+\beta PE{T}^2+\beta WeekPrecip+\beta SPEI+\beta Precip:AR+\beta SPEI:Precip+\beta SPEI:AR:Precip+\varepsilon .}$ We selected these variables to model hydroregimes because of typically strong ephemeral wetland responses to regional climate forcing (Brooks, 2004; Greenberg et al., 2015; Lee et al., 2015). We defined water inputs though direct precipitation effects and parameters that included precipitation within an interaction. The interaction between precipitation and the autoregressive term is useful for defining stage-dependent inputs of rainfall. We included linear and quadratic terms for PET to define short-term (i.e., daily) drawdown of water caused by evapotranspiration because this relationship was nonlinear at low PET values. Lastly, we included SPEI and parameters that include SPEI within an interaction to define wetland response to long-term drought conditions.

For each of the 35 well-specific models, we assigned vague priors following a normal distribution (mean = 0; variance = 100) to all fixed effects, with the error term assumed to follow a normal distribution with mean = 0 and the prior of variance to follow a uniform distribution (min = 0, max = 100). We fit each model using Markov chain Monte Carlo (MCMC) simulations, generating three chains, each with 400,000 total iterations and a thinning rate of 100 (Kéry & Royle, 2016). We used an adaptation phase of 1,000 and discarded 300,000 burn-in iterations, which retained 2,000 iterations for each chain (6,000 total samples) to estimate posterior distributions. We examined traceplots of parameters for adequate mixing among chains and the Gelman–Rubin’s $\hat{R}$ statistic to evaluate model convergence (Gelman, 2004). We assessed model predictive ability using a posterior predictive check based on the Bayesian P-value (Kéry, 2010; Link & Barker, 2010). We evaluated parameter significance based on the overlap of 95% highest posterior density with zero.

Model validation

For each well-specific model, we used 75% of the available data to train the model and 25% for testing. Because models were autoregressive, rather than using random data points, the training dataset included consecutive data starting at a random position within the first 25% of the data and consisted of 75% of these data while the testing dataset consisted of the remaining 25% of the data. To determine how well each model predicted the measured water level of our testing and training datasets, we calculated the normalized root mean squared error (NRMSE) for each of the 6,000 iterations, which is the root mean squared error divided by the range of measured water levels within each basin to account for differences in scale among basins. To determine how well each model predicted flooded conditions, we calculated the proportion of days in our testing and training dataset, for each iteration, that were predicted correctly to either have standing water or not.

Predicting hydrologic regime from future climate data

We obtained downscaled climate data for the years 2025–2100 using Localized Constructed Analogs (LOCA; Pierce, Cayan & Thrasher, 2014), which have a 6 × 6 km resolution. For each well, we obtained daily estimates of precipitation (mm), minimum temperature (°C), and maximum temperature (°C) for each of three global circulation models (GCM): Hadley Centre Global Environment Model 2 Earth Systems (HadGEM2-ES), Hadley Centre Global Environment Model 2 Carbon Cycle (HadGEM2-CC), and the Community Climate System model version 4 (CCSM4). For each model, we used two different projections for representative concentration pathways (RCP): 4.5 assumes a peak in carbon emissions in 2040 (Thomson et al., 2011), and 8.5 assumes emissions will continue to increase throughout the 21st century (Riahi et al., 2011). Therefore, for each model, we were able to assess future wetland hydrologic regimes under six different climate scenarios. For each of the six scenarios, we used the daily temperature and precipitation variables to calculate the same climate metrics (i.e., PET, SPEI, sum of rain over the last seven days) as the current data and predicted daily water level for each well using model parameter estimates from their respective posterior distributions.

To predict future hydrologic regimes and associated habitat suitability for flatwoods salamanders, we combined the LOCA climate data and the parameter posterior distributions from our statistical models. For each breeding season from 2025–2100, we determined the length (number of days) of the maximum hydroperiod within a given season, as well as the day in which each wetland became flooded. We defined the flatwoods salamander breeding season as occurring between 1 November and 31 May of the following year (Palis, 1996). Additionally, within a given posterior draw, we determined the number of consecutive breeding seasons that did not contain a hydroperiod of at least 15 weeks (105 days). While metamorphosis in flatwoods salamanders has occurred in as little as 11 weeks (77 days) and may take as long as 18 weeks (126 days; Palis, 1995), 15 weeks is a conservative estimate for the typical length of wetland filling required to allow larval development through metamorphosis (Brooks et al., 2020; Haas, 2010–2020, unpublished data) . Lastly, for each of the maximum hydroperiod estimates, we determined the date the wetland filled as the number of days since 1 November for the respective breeding season.

All analyses were performed in program R (version 4.1.1; R Core Team, 2021). We used the SPEI (Baguería & Vicente-Serrano, 2017) and Evapotranspiration (Guo, Westra & Peterson, 2020) packages for calculating SPEI and PET respectively, and the jagsUI package (Kellner, 2021) to call JAGS (Plummer, 2003), from Program R for MCMC analyses.

Model fitting

Our dataset consisted of 72,266 daily observations of wetland water level across 35 monitoring wells (range: 795–1,753 training observations per well; Table S1). These data covered a range of hydrologic conditions, both among wetlands (e.g., mean across wetland hydroperiod during the flatwoods salamander breeding season ranged from 37–142 days) and across years (e.g., mean annual breeding season hydroperiod across all wetlands ranged from 31–168 days). All MCMC chains showed good mixing, and $\hat{R}$ values were between 1.000 and 1.003, indicating model convergence. The posterior predictive check indicated that each model fit the data well (Bayesian P-value range: 0.49–0.51).

Our models generally showed good predictability of water level data, with the median training NRMSE ranging from 0.01 to 0.04 and the testing NRMSE ranging from 0.13 to 0.35 (Figs. S1A and S1B). Variance in NRMSE was always <0.01 for training and testing datasets and generally showed a decreasing trend with the number of training or testing datapoints. The median proportion of days correctly predicted to either have or not have surface water ranged from 0.94–0.99 for training datasets and from 0.56–0.93 for testing datasets (Figs. S1C and S1E). Averaged across all model iterations, the percentage of testing days predicted to be dry when the wetland was flooded ranged from 1.2–39% per wetland, while the percentage of testing days predicted to be flooded when the wetland was dry ranged from 0.3–29.4% per wetland. Some wetlands were more likely to misclassify wet or dry days, but others were equally likely to misclassify either condition. The variance in proportion of days correctly predicted to have surface water present ranged from <0.001–0.01 and from 0.12–2.40 for training and testing datasets, respectively (Fig. S1D and S1F).

Model parameters

The mean value for the first order autoregressive coefficient was positive and slightly less than one for all wetland basins. We found that daily precipitation had the largest magnitude effect, when compared to all other parameters, on wetland water level (Fig. 1; Table S2). The amount of precipitation over the previous seven days had a smaller positive effect than daily precipitation on water level. The model intercept (α), the SPEI, the quadratic term for PET, and the interactions between SPEI and daily rainfall and daily rainfall, autoregressive term, and SPEI had both negative and positive effects depending on the model. Finally, the effects of PET and the interaction between precipitation and the autoregressive term tended to have negative effects on daily water level (Fig. 1; Table S2).

Future predictions

Future projections revealed suitable hydroperiods for flatwoods salamander larval development (i.e., at least 15 consecutive weeks or 105 consecutive days of surface water during the breeding season) in at least some years across all climate scenarios (Figs. 2 and 3). When comparing the various climate scenarios, the CCSM4 model predictions typically had fewer years with suitable hydroperiods compared to the HadGEM2-CC or HadGEM2-ES model (Figs. 2 and 3). For date of wetland filling, wetland basins generally showed a tendency towards earlier fill dates (Fig. S2). For both hydroperiod and date of wetland filling, there was considerable variation among years and iterations for all wetland basins (Figures 2–3 and Fig. S2). Lastly, across all wetlands and climate scenarios examined in this study, we found substantial variability in their general suitability for flatwoods salamander reproduction (Fig. 4). We found that wetlands experiencing longer hydroperiods, on average, were less likely to have model iterations of at least five consecutive years of reproductive failure (i.e., five consecutive years with a hydroperiod less than 105 days). This metric is thought to be important for flatwoods salamander population persistence (Palis, Aresco & Kilpatrick, 2006). Some wetlands did have a high probability (>75%) of extended recruitment failures or had a low number of years with a suitable hydroperiod, indicating they are likely unsuitable for flatwoods salamanders. Overall, most wetland basins had suitable hydroperiods for flatwoods salamanders in at least 50% of future years (Figs. 2 and 4).