2.1. Study area

Our study area is coastal wetlands in the entire NGOM ranging from the western coastline of Florida in the east to the Texas coastline in the west. The majority of the NGOM is humid warm temperate/subtropical; although the southern tip of Florida, including the Florida Keys, has a tropical climate. Mangroves dominate along the southern coast of Florida (Mendelssohn et al. Citation2017) with Rhizophora mangle (red mangrove) dominating shorelines and creek banks, while monospecific stands of Avicennia germinans (black mangrove) or mixed stands of both species and/or Laguncularia racemosa (white mangrove) are in the landward areas. North of Tampa Bay, salt marsh gradually replaces mangrove (Osland et al. Citation2013) but due to its cold-tolerant trait (Markley, McMillan, and Jr Thompson Citation1982), only A. germinans is scattered along the northernmost region of the NGOM, including the Louisiana and Texas coasts (Mendelssohn et al. Citation2017).

The NGOM region contains the largest areas of salt marshes in North America with 55% of United States’ total (Mendelssohn et al. Citation2017). Salt marshes in western Florida, Alabama, and Mississippi are primarily irregularly flooded marshes dominated by Juncus roemerianus (black needlerush) with Spartina alterniflora (smooth cordgrass) at the seaward of the Juncus zone. Spartina patens (salt meadow cordgrass) and Distichlis spicata (salt grass) are at the higher elevations landward of the Juncus zone. Other common salt marsh species include Spartina cynosuroides(big cordgrass), Salicornia sp., Schoenoplectus americanus (chairmaker’s bulrush), and Schoenoplectus robustus (sturdy bulrush). Spartina alterniflora, A. germinans, and J. roemerianus dominate the deltaic marshes of Louisiana with subdominant species of S. patens and D. spicata. In the higher elevation of the Chenier Plain of southwestern Louisiana and southeastern Texas, the high marsh species, S. patens and D. spicata dominate brackish marshes. Salt marshes along the central Texas (Galveston Bay and adjacent bays) comprise mainly S. alterniflora, then Batis maritima (Saltwort), Salicornia spp., and J. roemerianus. The larger areas of brackish marshes at the higher elevations are dominated by D. spicata, Salicornia bigelovii (dwarf saltwort), Sarcocornia perennis (chickenclaws), Monanthochloe littoralis (shoregrass), and B. maritima with occurrence of S. patens and Spartina spartinae (gulf cordgrass). The marshes in the southern Texas transition to succulent-dominated hypersaline marshes, including B. maritima, Borrichia frutescens (sea ox-eye), Suaeda martima (herbaceous seepweed), and Sesuvium portulacastrum (sea purslane).

2.2. Data

We used five main data sources to develop the models. First, the USGS Hydrological Unit Codes 4 (HUC4) watershed boundaries (https://www.usgs.gov/core-science-systems/ngp/national-hydrography/watershed-boundary-dataset?qt-science_support_page_related_con=4#qt-science_support_page_related_con) to represent the pooled upland inputs (hydrology, sediments, and nutrients). It was used as a random factor to represent one broad spatial scale in the multi-level models (Figure 1). Second, the USGS streamflow data (ft³/sec) (https://waterdata.usgs.gov/nwis/sw), the only hydro-related data available at all the watersheds in the NGOM, used to represent watershed characteristics and the input to estuaries from uplands. Third, the USGS THK99 dataset (https://pubs.usgs.gov/dds/dds68/htmldocs/data.htm) for the geophysical covariates in the models including long-term average of RSLR rate, tidal range, mean wave height, and coastal slope (Figure 2). Fourth, Landsat-5 satellite images (https://earthexplorer.usgs.gov/) to calculate the normalized difference vegetation index (NDVI) to indicate vegetation productivity. Fifth, the NOAA Coastal Change Analysis Program (CCAP) data to 1) (https://coast.noaa.gov/digitalcoast/data/ccapregional.html) derive areal coastal wetland loss as the response variable in the models, 2) calculate percent developed areas including medium-high developed areas and cultivated lands, and 3) calculate percent forest areas including evergreen, deciduous, and mixed forests averaged for the years of 1996 and 2005 before Hurricane Katrina. Percent developed and forest areas represent watershed-scale anthropogenic influence.

Analyzing spatial variability of drivers of coastal wetland loss in the northern Gulf of Mexico using Bayesian multi-level models

All authors

Tyler Hardy , Wei Wu & Mark S. Peterson

https://doi.org/10.1080/15481603.2021.1943212

Published online:

01 September 2021

Figure 2. Boxplots of covariates and response variable for each watershed. The thick horizontal line shows the median, the wide box shows the upper and lower quartiles (25 and 75%), and the dotted “whiskers” show the distribution from the quartiles to the 1.5 interquartile range. RSLR denotes relative sea-level rise, and NDVI denotes normalized difference vegetation index from Landsat-5 imagery. Note: Landsat-5 images have seven spectral bands: band 1 = blue; band 2 = green; band 3 = red; band 4 = near-infrared (NIR); band 5 = shortwave infrared (SWIR 1); band 6 = thermal; band 7 = SWIR 2 with a revisit time of 16 days and a spatial resolution of 30 meters (m) except the thermal band with a spatial resolution of 120 m

Display full size

THK99 is a spatially-referenced dataset consisting of long-term (50–100 years) averages of RSLR rate, tidal range, wave height, and coastal slope at ~5 km intervals along the outermost shorelines of the contiguous United States and Alaska. The purpose of the data is to evaluate coastal vulnerability to relative RSLR (https://pubs.usgs.gov/dds/dds68/data/gulf/gulf.htm) and it has been applied to analyze coastal shoreline change (Thieler and and Hammar-Klose Citation1999; Gutierrez, Plant, and Thieler Citation2011a). The same drivers affecting shoreline change likely affect coastal wetland erosion. Furthermore, the shorelines the THK99 provides overlap with the boundaries of coastal wetlands, and therefore the data is appropriate for the analysis of vulnerability of coastal wetlands to RSLR. One limitation is that the wave data represents open estuarine conditions and may not accurately represent the wave height at coastal wetlands. However, this is the most consistent and best available wave data at the required region-specific scale and it has been used to classify vulnerability of coastal wetlands (Gutierrez, Plant, and Thieler Citation2011b).

We generated a point centroid for each 5-km polyline segment from the NGOM mainland shorelines (excluding barrier island shorelines) in the THK99 (n =1,184 segments), and created a circular buffer with a radius of 2.5 km around each centroid. We chose 2.5 km because it was half of the average distance between point centroids, which maximized spatial coverage while minimizing overlaps of buffered areas. We then used the buffer areas as the finest scale for modeling, described as “sites” in this paper. We extracted the geophysical properties at the centroid of each buffer area to represent the properties for that site area. If the buffer areas crossed two watersheds, we did not include them in the final analysis.

Watershed boundary represents a broader spatial scale than the site scale as each watershed contains multiple sites. We used 14 HUC4 watershed boundaries/hydrological regimes from the National Watershed Boundary Dataset (Figure 1) as a master variable to represent the pooled upland inputs, including freshwater, sediment, and nutrients, etc. We assumed the impacts of biogeophyscial variables on coastal wetland loss were more similar within the same watersheds than in different watersheds as the wetland sites in the same watersheds received similar upland inputs. We used the HUC4 to ensure there are enough sites within each watershed in order to obtain robust results on the impact of local drivers on coastal wetlands. We further evaluated whether inclusion of streamflow and land use/land cover at the watershed scales improved model predictions or not if the impact of local drivers differed by watershed.

The biological covariate used was NDVI, which is the normalized difference between the reflectance in the near infrared band and the red band. We downloaded the Landsat-5 images and calculated NDVI values using the Google Earth Engine (GEE) platform. The temporal coverage of the NDVI data was June through August of 1996, the season for peak vegetation productivity in the NGOM (Mo, Kearney, and Turner Citation2019; Mo, Kearney, and Momen Citation2017). We created a cloud-free image collection in the GEE for the entire temporal and spatial coverage and used the median NDVI of each 30 m pixel from June to August to generate the final composite image for NDVI. We extracted the NDVI values for only coastal wetland areas using GEE and then calculated the average of NDVI of the coastal wetlands for each study site. Larger NDVI values generally reflect higher aboveground green biomass (Couvillion and Beck Citation2013; Payero, Neale, and Wright Citation2004; Wen, Yang, and Saintilan Citation2012; Hardisky, Gross Gross, and Klemas Citation1986). We assumed healthy vegetation in 1996 helped maintain coastal wetlands for the next decade based on the observation that vegetation cover in prior year was the best single predictor of subsequent wetland loss in coastal Louisiana (Schoolmaster et al. Citation2018).

We derived the response variable, indicating areal wetland loss (ha) between 1996 and 2005, from NOAA’s CCAP data for the years of 1996 and 2005 (before Hurricane Katrina). CCAP land use/land cover data has a spatial resolution of 30 m, an appropriate resolution for a broad scale like the entire NGOM. It is likely this data may not capture coastal wetland loss finer than this resolution. We calculated areal wetland loss as the sum of cells, converted from any estuarine or palustrine wetlands to open water, at each site multiplied by the area of each cell. This estimate excluded the buffer areas without wetland loss that mainly represent areas with restoration activities or very limited coastal wetland areas. In total, we had 598 sites with wetland loss (Figure 1), and Figure 2 illustrates the five biogeophysical covariates and response variable by watershed.

2.3. Multi-level modeling in a Bayesian framework

Hierarchical Bayesian models allow the assimilation of data at multiple spatial scales through multi-level modeling and can quantify uncertainties (Clark Citation2005; Wu et al. Citation2012). They accommodate complexity by decomposing high-dimensional relationships into levels of conditional distribution within a consistent framework (Clark Citation2001; Wu, Clark, and Vose Citation2010). The framework includes a data level, a process model level, and a parameter level. In this framework, numerous latent variables and parameters that describe complex relationships are quantified in the form of posterior probability distributions that make uncertainty readily derived using credible intervals (Clark and Gelfand Citation2006). The Bayesian method usually applies Markov Chain Monte Carlo (MCMC) simulations to sample posterior distributions based on the product of likelihood and prior probability (Gilks, Richardson, and Spiegelhalter Citation1996).

Due to the multi-scale data (site and watershed scales) involved in analyzing coastal wetland loss, we used a multi-level modeling approach that is powerful when analyzing data based on nested designs (Bolker et al. Citation2009; Gelman and Hill Citation2007; Wu et al. Citation2018). We focused on studying the cross-scale interaction by examining how the effects of site-scale covariates on coastal wetland loss differed by watershed. The Bayesian inference accounts for uncertainties from data, processes, and parameters (Figure 3). Though some watersheds had only one site or a small numbers of sites, the variances could still be estimated due to the partial pooling with the multi-level models (Kéry Citation2010; Technow and Radu Totir Citation2015), but generally large. Partial pooling allows estimation of watershed-specific effects and the “borrowing” of information across watersheds. To represent the function of wetland loss (i.e., logarithm of wetland loss) at the watershed $i$ and site $j$ (W_ij), let ${\mathrm{W}}_{\mathrm{i}\mathrm{j}}\cdot$μ represent the mean of ${\mathrm{W}}_{\mathrm{i}\mathrm{j}}$, and ${\sigma }_{\hspace{0.17em}}^2$ represent the variance of wetland loss across sites. ${\mathrm{W}}_{\mathrm{i}\mathrm{j}}$ was modeled by assuming it was distributed as $(\sim )$a normal distribution $(N)$ (Equation 1): (1) $W_{ij}\sim N(W_{ij.}\mu ,{\sigma }^2).$(1)

Analyzing spatial variability of drivers of coastal wetland loss in the northern Gulf of Mexico using Bayesian multi-level models

All authors

Tyler Hardy , Wei Wu & Mark S. Peterson

https://doi.org/10.1080/15481603.2021.1943212

Published online:

01 September 2021

Figure 3. Structure of the Bayesian multi-level model used to identify the biogeophysical variables that affected wetland loss. Wetland loss was a function of geophysical properties and a vegetation index, the effects of which may $(covariatesk+1\dots m)$ or may not $(covariates1\dots k)$ vary in different watersheds. For the covariates that we assumed to affect wetland loss by watershed, we sampled their associated parameters for each watershed from the parameters at the scale of the entire study area (global scale). $m=totalnumberofcovariates$, $k=numberofcovariates$, the effects of which on wetland loss did not vary in different watersheds. $\beta s=interceptandparameters$ associated with normalized covariates $(C)$, $\alpha s=thehyper-parameters$ from which $\beta s$ were sampled, $W=logarithmofarealwetlandloss$, $i=watersheds$, and j = sites. Refer to the main text for formulas of different equations identified in this figure

Display full size

We modeled the mean wetland loss $(W_{ij.}\mu )$ using a linear function of the covariates C₁-C_m, with $k$ covariates having the same effect on wetland loss across the entire study area, and the remainder of the covariates $(m-k)$ having different effects on wetland loss in different watersheds. The specific functions representing the relationships between the potential drivers and coastal wetland loss are likely to vary by watershed, so we believe the linear representation of these relationships is necessary and appropriate at the scales we focused on for this study (sensu Braswell and Heffernan Citation2019). The intercept may vary $({\beta }_{i.,0})$ or not vary $(\beta .{.}_{,0})$ across watersheds. Assuming the intercept did not vary across the watersheds, we described the linear function using Equation 2: (2) $W_{ij}.\mu =f\left(\beta {..}_{,0},\beta {..}_{,1}\dots \beta {..}_{,k},{\beta }_{i.,\left(k+1\right)}\dots {\beta }_{i.,m}\right)={\beta }_{..,0}+\sum _{l=1}^k{\beta }_{..,l}{C}_{ij,l}+\sum _{l=k+1}^m{\beta }_{i.,l}{C}_{ij,l}$(2)

where βs represents the parameters for the linear model. If the intercept varied across watersheds, ${\beta }_{..,0}$ should be replaced with ${\beta }_{i.,0}$ in Equation 2. ${\beta }_{..,l}$ represents the coefficient for $C_{ij,l}$, for which the effects on wetland loss did not vary across the watersheds (1 ≤ l ≤ k). ${\beta }_{i.,l}$ represents the coefficient for $C_{ij,l}$, for which the effects on wetland loss varied across the watersheds (k + 1 ≤ l ≤ m). Therefore, we modeled wetland loss (W) for the R watersheds and S_i sites in each watershed (i) was modeled as in Equation 3: (3) $\begin{array}{r}p\left(\mathbf{W}|\beta ..{,}_0,\beta ..{,}_1\dots \beta ..{,}_k,{\beta }_{i\cdot ,}{ }_{\left(k+1\right)}\dots {\beta }_{i\cdot ,m}\right)\\ \propto \prod _{i=1}^R\prod _{j=1}^{S_i}N(W_{ij}|f\left({\beta }_{..,0},{\beta }_{..,1}\dots {\beta }_{..,k},{\beta }_{i.,\left(k+1\right)}\dots {\beta }_{i.,m}\right),{\sigma }^2).\end{array}$(3)

We sampled the coefficients that varied across the watersheds $({\beta }_{i.,\left(k+1\right)}\dots {\beta }_{i.,m})$ from the parameters at the coarser gulf-wide scale (or global scale in Figure 3) using normal distributions with means of αs and variances among watersheds of ${\eta }_{\hspace{0.17em}}^2$ (Equation 4): ${\beta }_{i.,\left(k+1\right)}\sim N\left({\alpha }_{..,\left(k+1\right)},{\eta }_{\left(k+1\right)}^2\right)$ ${\beta }_{i.,\left(k+2\right)}\sim N\left({\alpha }_{..,\left(k+2\right)},{\eta }_{\left(k+2\right)}^2\right)$

… (4) ${\beta }_{i.,m}\sim N\left({\alpha }_{..,m},\hspace{0.17em}{\eta }_m^2\right).$(4)

To complete the Bayesian model, we defined prior distributions for unknown parameters $(\alpha s,\beta s-{\beta }_{..,1}to{\beta }_{..,k},{\sigma }^2,and{\eta }^{2}s)$. We used conjugate prior for computation efficiency (Calder et al. Citation2003) and therefore the priors and posteriors had the same probability distribution forms. The priors for $\alpha s$ and $\beta s$ were normally distributed, and the priors for ${\sigma }^2$ and ${\eta }^{2}s$ followed the inverse gamma distributions. The prior distributions were flat and only weakly influenced the posteriors, which reflected the lack of knowledge of these parameters (Lambert et al. Citation2005) (see supplementary material for the specification of the priors).

By combining the parameters (priors), processes, and data models, we derived the joint posterior distribution of unknowns (Clark Citation2005) in Equation 5: (5) $p({\beta }_{\mathrm{..},0},{\beta }_{\mathrm{..},1}\dots {\beta }_{\mathrm{..},k},{\beta }_{i.,\left(k+1\right)}\dots {\beta }_{i.,m}{\alpha }_{\mathrm{..},\left(k+1\right)}\dots {\alpha }_{\mathrm{..},m},{\sigma }_{ }^2,{\eta }_{\left(k+1\right)}^2,\dots ,{\eta }_m^{2}W,{\text{C}}_{\text{1}}\dots C_{\text{m}})$(5) $\begin{array}{rl}& \propto p\left(\mathit{W}|f\left({\beta }_{..,0},{\beta }_{..,1},\dots ,{\beta }_{..,k},{\beta }_{i..,\left(k+1\right)},\dots ,{\beta }_{i..,m}\right),{\sigma }^2\right)\\ & \times \left(p\left({\beta }_{i.,\left(k+1\right)}|{\alpha }_{..,\left(k+1\right)},{\eta }_{\left(k+1\right)}^2\right)\times \dots \times p\left({\beta }_{i.,m}|{\alpha }_{..,m},{\eta }_m^2\right)\right)\end{array}$$\left.\begin{array}{rl}& \times p\left(\beta |..,0\right)\times \left(p\left({\beta }_{..,1}\right)\times \dots \times p\left({\beta }_{..,k}\right)\right)\times p\left({\sigma }^2\right)\\ & \times \left(p\left({\alpha }_{..,\left(k+1\right)}\right)\times \dots \times p\left({\alpha }_{..,m}\right)\right)\\ & \times \left(p\left({\eta }_{\left(k+1\right)}^2\right)\times \dots \times p\left({\eta }_m^2\right)\right)\end{array}\right\}\mathrm{P}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}\mathrm{s}$

The model in Equation 5 was parameterized using MCMC simulations (Gelfand and Smith Citation1990) in JAGS through the package “rjags” (Plummer Citation2016) in CRAN R (R Core Team Citation2015) (see supplementary material for the JAGS code). Prior to implementation of the multi-level models, we performed data transformations to improve the model computation efficiency. The covariates were normalized using classic standard score normalization (Zar Citation2010). Normalization allowed comparison of the covariates’ impact on wetland loss by examining the magnitudes of their coefficients directly because they were transformed to the same ranges (Hooten and Hobbs Citation2015). Furthermore, the response variable was log transformed to fit a normal distribution and account for nonlinearity and non-additive effect of the biogeophysical variables. We checked for multicollinearity using variance inflation factor (VIF). We only removed a covariate from further consideration if the VIF for it was greater than five (Zuur et al. Citation2009). None of the covariates showed VIF ˃5, so we did not remove any of the covariates.

We simulated three MCMC chains using three different sets of initial values (based on vague priors) for the parameters. The three chains were used to determine the number of iterations of MCMC before they converged for the model (20,000 iterations) (Cowles and Carlin Citation1996). We then discarded the pre-convergence, burn-in iterations, and ran the chains additional 300,000 iterations thinning every 10 steps to reduce within-chain autocorrelation. The generated posteriors helped identify the covariates that had effects on wetland loss different from zero, evaluate relative magnitude of the effect of each covariate on wetland loss, and determine whether a covariate had different effects on wetland loss in different watersheds.

2.4. Model comparison

We compared the 484 models, based on all combinations of covariates and whether they interacted with watersheds or not, to identify the best model(s) using both deviance information criterion (DIC, Eq. 6) and posterior predictive loss (PPL, Eq. 7) functions. The lower the DIC or PPL, the better the model predicts (Hooten and Hobbs Citation2015). Although PPL functions measure the error in model predictions as loss (Gelfand and Ghosh Citation1998; Ibáñez et al. Citation2009), DIC generally outperforms PPL in comparing models with large sample size (n >≈ 1000) (Daniels, Chatterjee, and Wang Citation2012). However, DIC only allows for comparison of models with the same level of structures (Hooten and Hobbs Citation2015). Furthermore, models are considered comparable if the difference of their DICs is less than 2 (Burnham, Anderson, and Huyvaert Citation2011; Burnham and Anderson Citation2004). Therefore, we used DIC to compare models with the same hierarchies (i.e., comparing models considering watersheds’ effect but with different covariates) and PPL to compare the models with different structures (i.e., comparing models considering watershed’s effect to those without). We identified the best model with the lowest PPL and DIC, while considering models as the best candidate models if their DICs are not larger than the DIC of the best model by more than 2 (Hobbs and Hooten Citation2016). (6) $DIC=-2\mathrm{l}\mathrm{o}\mathrm{g}([\mathit{W}|E\left(\mathit{\theta }|\mathit{W}\right])+2p_D$(6)

where $p_D$ represents an effective number of parameters, θ represents the parameters in the model, E represents expected values, and [] represents probability. (7) $PPL=\sum _{i=1}^n{\left(y_i-E\left(W_i^{new}|\mathit{W}\right)\right)}^2+\sum _{i=1}^{n}Var\left(W_i^{new}|\mathit{W}\right)$(7)

where $W_i^{new}$ represents a new data point (wetland loss) simulated from the model, and Var represents variance.

If the best model(s) showed the impact of the covariates or intercept varied by watershed, we further added mean daily streamflow (ft³/sec), and/or land use/land cover as covariates at the watershed scale and investigated whether the inclusion of them improved the model predictions. In the models that accounted for these watershed-scale covariates, we modeled the mean of ${\beta }_{i.,0}\left({\beta }_{i.,0.}\mu \right)$as a function of streamflow or streamflow and land use/land cover (Equations 8–10). There was no multicollinarity problem when we integrated streamflow and land use/land cover as the highest correlation was −0.37 between percent forest areas and percent developed areas. We compared these different models to the best model(s) based on DIC only as they had the same hierarchies. (8) ${\beta }_{i.,0}.\mu =k_0+k_{1}SF,$(8) (9) ${\beta }_{i.,0}.\mu =k_0+k_{1}SF+k_{2}DV+k_{3}FR,$(9) (10) ${\beta }_{i.,0}\sim N\left({\beta }_{i.,0}.\mu ,{\eta }_0^2\right)$(10)

where k₀-k₃ represent the parameters, SF represents streamflow, DV represents % developed areas, FR represents % forest areas at the watershed scales, and ${\eta }_0^2$ represents variance of ${\beta }_{i.,0}$.

We examined the 95% credible intervals (CIs) of the posteriors for the coefficients to determine whether the covariates had effect on wetland loss different from zero (not overlapping 0) with the sign of the CI for a given covariate indicating a positive (CI is positive without overlapping zero) or negative (CI is negative without overlapping 0) effect on wetland loss.

We considered Bayesian p-values (Equation 11) as a method to examine the model fitting in moments and quantiles. Bayesian p-values describe the proportion of times the statistic of interest generated from the model is equal to or more extreme than the statistic from the data; thus, they differ from traditional p-values. The value of 0.5 indicates a good fit of the model to the data while the values close to 0 or 1 indicate a lack of model fit that cannot be attributed to chance. (11) $P=\mathrm{p}\left[T\left(y_{rep},\theta \right)\ge T\left(y,\theta \right)|y\right]or=\mathrm{p}[T\left(y_{rep},\theta \right)\le T\left(y,\theta \right)|y]$(11)

where $P$ denotes p-value, $T\left(y,\theta \right)$ is a test-statistic of interest (i.e., mean, standard deviation, quantiles) calculated from observed data (y), while $T\left(y_{rep},\theta \right)$ is a corresponding statistic from the new “data” $\left(y_{rep}\right)$ generated from the posterior predictive distribution. These estimates allow for the calculation of p-values for various statistics (Meng Citation1994). Table S1 provides all of the structures, DICs, PPLs, and Bayesian p-values for the 484 models.

Multi-scale considerations

In the limited research of coastal wetland dynamics that account for multiple spatial scales, Braswell and Heffernan (Citation2019) focused on wetland distribution at the estuary and watershed scales instead of wetland loss studied here. It did not address or quantify the variance of the impact of the biophysical factors on coastal wetlands either. Their results are, however, useful as they indicated that characteristics of estuaries (depth, shoreline complexity, land use, tidal range, and latitude) and their watersheds (discharge) best predicted wetland extent at broad spatial scales (~>10² km²). In their model, bimodal distribution of wetland extent occurs at finer scales (~10⁰ to 10² km²) that correspond to the theoretical scale of wave erosion feedbacks. These results suggest that local feedbacks are the mechanism that controls marsh dynamics, whereas coastal and watershed characteristics determine macroscale wetland distributions through their effects on these local processes.

The interactions between the two spatial scales selected in our final models pinpoint the importance of accounting for watershed boundaries in evaluating coastal wetland loss as upland and marine processes influence coastal wetlands in transitional environments. The coastal sites within the same watershed tend to receive similar upland inputs like decreased freshwater, pollutant discharges, land use/land cover associated changes, and changing climate responses compared to sites in the other watersheds (Hovenga et al. Citation2016; Alber Citation2002; Lin and Shen Citation2018; Huang et al. Citation2016; Li et al. Citation2018). Therefore, the watershed boundary in our model coarsely reflects the conditions of upland inputs. The changes in upland inputs affect salinity, nutrients, and temperature dynamics in the systems (Kennedy and Barbier Citation2016; Elsey-Quirk et al. Citation2019) and further influence the health of coastal wetlands in the short and long terms. Importantly, upland watershed modifications reduce or eliminate river-borne sediments which are the main contributor of vertical accretion, formation, and maintenance of coastal wetlands in the NGOM or elsewhere (Karegar, Dixon, and Malservisi Citation2015; Tweel and Turner Citation2012; Weston Citation2014; Blum and Roberts Citation2009). In fact, the fluvial sediment input to the NGOM substantially decreased from 1950 to 2010 due to reduction in agricultural soil erosion (Syvitski et al. Citation2005), damming of rivers (Vörösmarty et al. Citation2003), and water-use in the upper watershed (Meade and Moody Citation2010), which in turn causes extensive loss of coastal wetlands in the NGOM. Because sediment input data is limited and inconsistent across our large study area, including streamflow and land use/land cover as covariates at the watershed scale likely captures some of the sediment load and nutrient inputs that would represent upland inputs and thus may influence our model results. However, when we accounted for streamflow and land use/land cover in the linear models at the watershed scale, the model predictions did not improve. This is likely due to these two covariates not correlating well with sediment load due to hydrological alterations spatially and temporally on the rivers. The models should also consider nonlinear functions at the watershed scales.

Drivers for coastal wetland loss

The coefficient for RSLR had both a positive median value for each watershed except for southern Florida (#1). RSLR also had the largest median effect on wetland loss among all the covariates in the majority of watersheds (Figure 4). These results are supported by studies that indicated an increase of wetland loss due to increasing sea level (Kirwan and Patrick Megonigal Citation2013; Linhoss et al. Citation2015; Morris et al. Citation2002; Simas, Nunes, and Ferreira Citation2001; Wu et al. Citation2015; Wu, Biber, and Bethel Citation2017). However, large spatial variability existed in the impact of RSLR on areal wetland loss across the NGOM. The north-central and north-eastern Gulf of Mexico, including Choctawhatchee-Escambia watershed in western Florida panhandle (#6), Mobile-Tombigbee watershed in Alabama (#7), Pascagoula watershed on MS coast (#8), MS Delta watershed in Louisiana (#9), and LA coastal (or western Louisiana or Chenier Plain) watershed (#10) from east to west, showed highest vulnerability to RSLR based on the posteriors of the coefficient. The patterns we identified are supported by research that indicated high RSLR caused increases in coastal wetland inundation (Alizad et al. Citation2016), non-linear expansion of tidal creeks (Hagen et al. Citation2012), and coastal wetland loss (Geselbracht et al. Citation2015; Glick et al. Citation2013; Kidwell et al. Citation2017; Martin and White Citation2000; Passeri et al. Citation2016; Walsh Citation2007).

We expected the RSLR coefficients to be positive and largest for both the MS Delta watershed (#9) and LA coastal (or Chenier Plain) watershed (#10) due to the extensive literature on the subject. For example, high rates of wetland loss was due to high RSLR from sediment compaction and subsidence in Louisiana (Georgiou, FitzGerald, and Stone Citation2005; Karegar, Dixon, and Malservisi Citation2015; Mallman and Zolback Citation2007; Törnqvist et al. Citation2008; Yuill, Lavoie, and Reed Citation2009; Couvillion et al. Citation2017; Bourne Citation2000). In fact, the medians of the coefficients for RSLR in both Louisiana’s watersheds were positive and RSLR showed positive and marginally positive effects in the LA coastal (or Chenier Plain) watershed (#10) and MS Delta (#9) watersheds, respectively. The RSLR’s median effect on wetland loss in the LA coastal (or Chenier Plain) watershed was larger than in all other watersheds. This result was consistent with the previous research that showed the LA coastal (or Chenier Plain) watershed was more vulnerable to RSLR compared to the MS Delta due to differences in geomorphological features and processes like proximity to riverine sediment inputs, connectivity to the Gulf of Mexico, impact of Chenier ridges and impoundments (Jankowski et al. Citation2017).

The coefficient for wave height (Figure 4) was positive only in the MS Delta watershed (#9), indicating that wave height had a positive effect on wetland loss that was comparable to that of RSLR. The Mississippi Delta is characterized by high rates of subsidence partially due to compaction caused by decomposition of carbon-rich Holocene era sediment deposits (Törnqvist et al. Citation2008), large land loss (Hartig, Kolker, and Gornitz Citation2002; Kolker, Allison, and Hameed Citation2011), and/or pores left from oil extraction (Kolker, Allison, and Hameed Citation2011; Mallman and Zolback Citation2007). The sediments of the MS Delta in Louisiana are also less compacted than sediments in most of the rest of the NGOM due to larger total porosity (Lisle and Comer Citation2011), and thus are more likely to be eroded by wave action (McLoughlin Citation2010; Fagherazzi et al. Citation2013; Allen Citation1989). Thus, increased wave height in the MS Delta watershed would lead to increased wetland loss as derived by our model, irrespective of RSLR.

In contrast, the negative effect of wave height on coastal wetland loss in the southern Florida watershed (#1) is unexpected but is likely purely correlational due to low wave height and high wetland loss in this watershed. It may be better interpreted as that coastal wetland loss was large even under low wave height (Figure 2). Mangroves in the southern Florida watershed instead of salt marshes may explain low wave height as they reduce wave energy and wave height because of mechanical friction from roots and trunks (Parvathy and Bhaskaran Citation2017). Mangroves loss is large in this watershed, likely due to complex anthropogenic influences beyond land use/land cover change (Goldberg et al. Citation2020).

It was expected that coastal slope and tidal range would lower coastal wetland loss as steeper coastal slopes act as a physical barrier to inundation and larger tidal range could deposit more ocean-borne sediment onto coastal wetlands (Kirwan and Guntenspergen Citation2010). The general negative medians of their coefficients for each watershed in our study support this assertion (Figure 4). However, a steeper slope does not facilitate upslope migration of salt marshes and thus may contribute to wetland loss. We did not consider this process in our model, but the wetland change may show up in the response variable. The dual effects of coastal slope may explain why it was not selected in the best model (Model 58).

NDVI did not affect coastal wetland loss shown in the models as it was either not included in the best model candidates or the posterior of its coefficient overlapped zero. This may be because NDVI may not adequately represent vegetation productivity. Although NDVI has been used in the NGOM to predict above-ground vegetation productivity and biomass (Couvillion and Beck Citation2013; Mishra et al. Citation2012; Mo, Kearney, and Momen Citation2017; Mo, Kearney, and Turner Citation2019), an evaluation of multiple remote sensing derived vegetation indices indicated a large spatial variability in the index that best represents peak green above-ground biomass (Wu et al. Citation2018). Some studies show wetland health is rarely correlated with above-ground biomass (Nyman et al. Citation1994; Turner et al. Citation2004; Kirwan and Guntenspergen Citation2012) as below-ground biomass is key to stability of coastal wetlands (Tahsin, Medeiros, and Singh Citation2016; Byrd et al. Citation2018; Turner et al. Citation2004; Stagg et al. Citation2016). The characteristics of below-ground biomass cannot be directly related to surface reflectance detected by remote sensors. Inundation levels in coastal wetlands can also affect NDVI and NDVI-based estimates of biomass (Klemas Citation2013; Kearney et al. Citation2009). Finally, our choice of the 1996 period (June – August) for NDVI, the beginning of the decadal change, may not be ideal to account for the long-term effect of vegetation status. Though vegetation cover in the prior year was the best single predictor of subsequent coastal wetland loss in coastal Louisiana (Schoolmaster et al. Citation2018), this does not necessarily extrapolate to a long-term and larger spatial extent investigation as done in our study. NDVI in the middle of the decade or the mean decadal NDVI may be worth exploring further.

Missing drivers

We did not select all the potential drivers in our study as illustrated by the Bayesian p-value (based on coefficient of variance) of 0.28 for the best model. We chose to focus on potentially important biogeophysical drivers that are readily available spatially. Some fine-scale covariates (i.e., site scale) that can affect coastal wetland vulnerability to RSLR are salinity, river/tide ratio (Braswell and Heffernan Citation2019), interior marsh inundation (Day et al. Citation2000), wind driven wave height (Marani et al. Citation2011), vegetation type, and associated flow attenuation (Rodríguez et al. Citation2017). In addition, storms, sediment deprivation, oil and gas extraction and infrastructure, navigation infrastructure, and saltwater intrusion altered hydrology and likely coastal wetlands (Penland et al. Citation2001; Couvillion et al. Citation2017). For example, river borne-sediment deprivation in Louisiana as indicated by a 50% decline in sediment supply since 1950s due to upland hydrological changes (Reyes et al. Citation2000; Blum and Roberts Citation2009), progressive soil waterlogging due to reduced sediment input, and interaction with salinity, produce deleterious effects on vegetation (Reyes et al. Citation2000; Turner Citation1990).

Anthropogenic factors were not explicitly accounted for in our models of wetland loss except land use/land cover (Lin and Shen Citation2018; Li et al. Citation2018), though some impacts were likely pooled into watershed boundaries, such as pollutant discharge (Lin and Shen Citation2018), damming, dredging and navigation activities (Turner Citation1997; Li et al. Citation2018), and population density (Entwistle, Mora, and Knight Citation2018). However, some human activities occur at the estuary scale (i.e., between site and watershed scales) and the addition of the estuary scale would likely improve the models further. We also did not include climate variables and soil properties directly in our models. However, climate variables directly affect vegetation growth, so the vegetation index would likely capture some of the variability of climate variables. Furthermore, RSLR may affect wave height or tidal range (Passeri et al. Citation2016), while their impacts to each other were not included in the current study. Finally, we used HUC4 as a balanced approach between watershed size and number of sites in each watershed, and we believe the choice of scales in our study is appropriate. However, the selection of watershed size can be more refined than HUC4 if a sufficient number of sites are available in each watershed. More covariates like sediment availability that are relevant to the study design are important to be included as they should improve the model predictions if they are available.

Implication for future research and resource management

The Bayesian multi-level modeling approach we developed provides a flexible framework that can be easily adapted to incorporate more data while updating the posteriors with uncertainty reduced. It has the flexibility of integrating different spatial-scale data (see Figure 3, such as adding a new spatial scale, or adding covariates at an existing spatial scale like indicators for anthropogenic activities), and can be applied more broadly in other regions with different environmental settings as long as the required data on wetland loss and potential drivers are available. This modeling approach is useful in identifying the biogeophysical factors affecting wetland loss and detecting hotspots of vulnerability (tier 1). It facilitates more in-depth analysis at the hotspots of vulnerable wetlands that would require urgent attention and timely collection of more fine-scale data to understand the processes that contribute to that vulnerability (tier 2). We found we could reduce research labor costs and resources by using the two-tier research approach but the modeling framework may be limited by the lack of data for some fine-scale drivers.

In addition to increased understanding of coastal wetland dynamics, quantifying the dominant drivers for coastal wetland loss at a particular watershed, as done in our study, has important implications for coastal wetland restoration. In watersheds where dominant drivers are difficult to identify, more comprehensive restoration strategies that address multiple drivers will be needed. If RSLR is the dominant driver, sediment enhancement through beneficial use of dredged materials and/or river diversion, or increases in hydrological connectivity may be effective (Ravens et al. Citation2009), like in MS coast (#8), western Louisiana or Chenier Plain (#10), and MS Delta (#9) watersheds. If wave height is the main driver, breakwater structures or living shorelines to reduce wave energy may be considered to mitigate wetland loss more effectively like in MS Delta (#9) (Davis et al. Citation2015). Finally, the real scenarios may require multiple strategies simultaneously to conserve or restore coastal wetlands effectively due to the complexity of coastal wetland dynamics. For the MS Delta (#9), both RSLR and wave height positively (or marginally positively) affected coastal wetland loss and thus effective restoration strategies would require the combination of freshwater diversion, living shorelines and/or breakwater infrastructure (Hardy Citation2018).