Climate drivers and human impacts shape 35-year trends of coastal wetland health and composition in an urban region

The future of coastal wetlands will depend on the combined effects of climate change and human impacts from urbanization and coastal management. Disentangling the effects of these factors is difficult, but satellite imagery archives provide a way to track biological and physical changes in wetlands over recent decades to reveal how coastal wetlands have been changing in response to climate and human drivers. In this study, we used Landsat to monitor the conditions of 32 coastal wetlands in southern California from 1984 to 2019 and identify environmental and human drivers of these trends. Wetland conditions were characterized by vegetation greenness, using the normalized difference vegetation index (NDVI)


INTRODUCTION
Coastal wetlands are some of the world's most valuable ecosystems because of their ability to store carbon, protect coastlines, and provide many other benefits to humans (Barbier, 2019;Barbier et al., 2011).Across the world, these ecosystems are threatened by climate change and human impacts, and it is unclear how coastal wetlands may respond to changing environmental conditions and to ongoing human influence (Kirwan & Megonigal, 2013).Understanding the relative importance of environmental drivers and human impacts on coastal wetlands, and how driver influence may vary over space and time, will be essential to adaptation efforts aimed at preserving these valuable ecosystems in the future (Erwin, 2009;Levin, 1992).
Coastal wetlands are sensitive to changes in climatic and environmental drivers (Osland et al., 2016).Rising seas, CO 2 concentrations, and temperatures, extreme storms and tidal events, and altered precipitation regimes are already impacting wetland ecosystems across the world (IPCC, 2022;Mckee et al., 2012), eliciting responses that can range from widespread habitat loss or degradation to shifts in wetland plant communities, as well as changes in plant productivity and phenology (Morris et al., 2002;Short et al., 2016).Response to environmental change can vary greatly among wetland species, sites, or systems due to biogeomorphic processes that balance sediment, fresh and salt water, and plant productivity (Janousek et al., 2016;Scavia et al., 2002).
Direct human impacts in coastal regions also significantly affect wetland ecosystems (Gedan et al., 2009;Newton et al., 2020).The urbanization of coastlines worldwide has led to wetland habitat loss and degradation, increased impervious surface cover, and shoreline hardening (Gittman et al., 2015).Hardened infrastructure can increase the risk of coastal squeeze (Borchert et al., 2018), habitat loss (Gittman et al., 2015), and cause changes to wetland plant communities (Watson et al., 2017).Additional stressors that often accompany urbanization include resource extraction, species introductions, and invasive species, hydrologic alteration, pollution, and eutrophication, and the impacts of urbanization are often highly contextual and vary geographically (Gedan et al., 2009).Coastal management and restoration programs have arisen in recent decades to counteract the negative impacts humans have historically had on wetlands (e.g., Callaway et al., 2011;Kennish, 1999;Stein et al., 2020;Zedler & Callaway, 1999); however, these efforts can also have direct, and sometimes unintended, consequences to plant communities, habitat distributions, and conditions within wetland and estuarine systems (Kennish, 2001).Sediment augmentation, for example, can bolster elevations in sinking marshes but has shown mixed effects to plant recovery and communities in wetlands of Louisiana, North Carolina, New York, and California, USA (McAtee et al., 2020).We can begin to uncover the comingling effects of human impacts on wetlands, in addition to environmental drivers, by tracking wetland health and restoration trajectories.
The coastline of southern California, in contrast, has unique physical features, climate, and a history of urbanization that has created a diverse network of over 80 fragmented coastal wetlands that range in size, setting, ecology, and typology (SCWRP, 2018).Wetlands remaining today are relatively small or discrete remnants comprising just 25% of areas that existed ca.1850 due to land use conversion (Grossinger et al., 2011;Stein et al., 2014Stein et al., , 2020) ) and shoreline hardening (Griggs & Patsch, 2019).Wetlands here experience a semiarid, Mediterranean climate marked by dry, hot summers and cool, wet winters, making them sensitive to changes in hydrology.The future of southern California's diverse wetlands is threatened by sea level rise (SCWRP, 2018;Thorne et al., 2016Thorne et al., , 2018)), and although up to 48% wetland habitat loss is predicted regionally, there is high uncertainty among wetland sites (Doughty et al., 2018).How individual wetlands have responded to both human and environmental stressors over the recent past can provide vital information to coordinated regional efforts managing, restoring, and adapting coastal wetlands into the future (e.g., SCWRP, 2018;Zedler, 1996).
In this study, we track the past conditions of 32 distinct wetland sites across the southern California region from 1984 to 2019 using the Landsat archive.Wetland conditions were measured as greenness (normalized difference vegetation index [NDVI]), an indicator of the biological and physical properties within wetlands, and as habitat composition, which we derived from the areal estimates of wetland (vegetated marsh and unvegetated mudflat) and subtidal habitats.Our aim is to (1) quantify long-term trends and variability in wetland cover and greenness over 35 years, and (2) identify the relative importance of environmental drivers to conditions in individual wetlands and for the region.We predict that the dominant environmental controls to past changes in wetland NDVI are related to precipitation, a major constraint in Mediterranean climates, and that changes to wetland habitat composition over this time period are driven by sea level rise.To understand human impacts on wetland conditions, we also aim to (3) compare the functional response (i.e., combined NDVI and habitat trends) to metrics of human impacts like urbanization and restoration activities for each site.This novel application of the Landsat archive in a region facing a combination of climate and human stressors can help uncover the drivers of future wetland response to environmental change and inform future management decisions.

Description of study region, wetland sites, and archetypes
Coastal wetlands in southern California are structurally and functionally diverse due to coastal setting and exposure (Jacobs et al., 2011), ranging in size and function from vernal pools occupying a few hundred meters to large estuaries over 250 ha (Stein et al., 2020).A wetland archetype typology was previously developed to classify wetlands of similar setting, structure, function, and plant community composition and to provide a framework to predict future response to SLR and aid ongoing regional restoration planning efforts (Doughty et al., 2018;SCWRP, 2018;Stein et al., 2020).Wetland "archetypes" used in this study include small creeks and lagoons (6), intermediate estuaries (6), large lagoons (7), large river valley estuaries (7), and fragmented river valley estuaries (6).We selected these 32 wetland study sites for a representative sample of archetypes that met the thresholds of size, shape, and habitat composition for Landsat analysis, and that have data available on management, restoration history, and environmental sources (Appendix S1: Table S1).
Our wetlands sites were manually delineated by the Southern California Coastal Water Research Project (SCCWRP) using 2016 National Wetland Inventory (NWI) data and aerial imagery.We used the 2016 SCCWRP wetland boundaries as the basis for this study, as these align with SLR-response modeling (Doughty et al., 2018) and adaptation planning and management guidelines (SCWRP, 2018).Boundaries used for each wetland contain contiguous, interconnected habitats of vegetated marsh, unvegetated mudflats and salt flats, and subtidal habitats (SCWRP, 2018).

Landsat time series of wetland greenness and habitat composition
We developed Landsat-based time series of wetland greenness and area for each of the 32 wetlands using Google Earth Engine (GEE), a cloud-based platform for storing and analyzing geospatial and remotely sensed data (Gorelick et al., 2017).Using GEE, we accessed Landsat Collection 1 Tier 1 Surface Reflectance scenes from 1984 to 2019 collected with Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI sensors.Landsat scenes contain pixels with 30-m ground spatial resolutions and ≤12-m root mean square error in geometric accuracy (Wulder et al., 2019).To improve continuity in the 35-year time series, we corrected for radiometric differences across sensors by harmonizing Landsat data using correction factors for surface reflectance (Flood, 2014;Roy et al., 2016;Vogelmann et al., 2016).We used the Landsat pixel quality assessment product ("Pixel_QA," formerly "CFmask"; Zhu & Woodcock, 2012) to remove pixels containing clouds and cloud shadows in each scene.We removed scenes from the time series if the total clear pixels remaining (including both water and land) was less than 75% of the given wetland area to avoid scenes with clouds and Landsat 7 ETM+ Scan Line Corrector (SLC) Failure (Andrefouet et al., 2003).
The resulting filtered collection of cloud-masked Landsat scenes was used to estimate three metrics of wetland conditions over time: the area of subtidal habitat, the area of wetland habitats (vegetated marsh and unvegetated mudflats), and the NDVI (Rouse et al., 1974) of wetland habitats.For area estimates, we developed additional masks from the Pixel_QA bitmask band that correspond to wetland habitats (bit 1: clear land) and subtidal habitats (bit 2: water).The area of each habitat was estimated as the sum of pixels in each mask multiplied by pixel area.We used this approach to consistently monitor site conditions and track changes over time in the relative distributions of wetland habitats.The CFmask has shown 89% accuracy in estimating marsh areas (Mo et al., 2019).
Only pixels identified as clear wetland habitats were used to calculate site NDVI as a metric of wetland greenness that is distinct from wetland area.Our previous work in the coastal wetlands of southern California region showed that NDVI was positively correlated with aboveground biomass (AGB) in simple vegetation canopies and AGB did not saturate with NDVI (Doughty et al., 2021;Doughty & Cavanaugh, 2019).Thus, we use NDVI as a proxy to quantify trends indicating biological and physical changes: greening reflecting increases in plant productivity, AGB, or vegetated cover in a wetland, and browning indicating habitat loss, conversion to mudflat or subtidal, or reductions in plant productivity or biomass.
Annual time series of wetland conditions, including wetland area, subtidal area, and NDVI in wetland habitats, were created for each of the 32 sites from 1984 to 2019 (Appendix S1: Figure S1).We performed additional quality control measures to inspect Landsat scenes with low pixel counts and low or outlying NDVI values to remove scenes with errors in cloud masking and filtering from the annual time series (Appendix S2).Landsat scene availability per site and year ranged from 0 to 65 images with a median of 15 images year −1 (Appendix S1: Figure S1).Time series filtering has been shown to reduce bias toward low NDVI and reduce variance in the data (Younes et al., 2021) and has no negative impacts on marsh phenology analyses (Buffington et al., 2018).We aggregated our time series by the mean on monthly and annual timescales to help average out the effects of tidal variations among images.We developed an interactive GEE application to allow users to view time series for the coastal wetlands across southern California (cheryldoughty.users.earthengine.app/view/so-cal-wetland-health).

Environmental and climatic datasets
To investigate the drivers of wetland conditions over time, we collected time series of environmental and climatic data from 1984 to 2019 that fall into seven categories: sea levels, wave conditions, stream discharge, precipitation, temperature, drought, and climate oscillations (Figure 1, Table 1; Appendix S2).Daily sea level data were taken for each site from the nearest available of five NOAA Stations in the region (Appendix S1: Figure S2).Data on daily stream discharge were available from USGS Stream Gauges for 11 of the 32 study sites (Appendix S1: Figure S2).
We used modeled wave data from the Coastal Data Information Program (CDIP) as time series of nearshore wave conditions at each site (O'Reilly et al., 2016), and sites were grouped into five wave climate subregions, the Gaviota Coast, Ventura Coast, Santa Monica (SM) Bay, San Pedro (SP) Bay, and San Diego (SD) Coast (Appendix S1: Figure S2), as these subregions reflect changes in regional topography (steep vs. shallow grades) and coastal orientation (southerly vs. westerly facing) that determine coastal wave exposure (Jacobs et al., 2011;SCWRP, 2018).Because the effects of waves play an important role in estuary mouth dynamics of intermittently opening and closing estuaries in this region (Clark & O'Connor, 2019), each wetland was categorized as either predominantly open (<40%), intermittently open/closed (40%-60%), or predominantly closed (>60%) based on the estimates of percent closure over time derived from the California Coastal Records Project (Doughty et al., 2018).
Monthly time series of temperature, precipitation, vapor pressure deficit (VPD), the self-calibrating Palmer drought severity index (PDSI), the standardized precipitation-evapotranspiration index (SPEI), and the standardized precipitation index (SPI) for each site were derived from gridded (4-km) PRISM data (Daly et al., 2008).We also considered climatic fluctuations using monthly data on the Oceanic Niño Index (ONI), the Pacific Decadal Oscillation (PDO) index, and the multivariate El Niño-Southern Oscillation (ENSO) index.
Environmental data specific to each wetland site were aggregated into annual time series and combined with the corresponding annual time series of wetland conditions for each site.All environmental and climatic variables were aggregated using the mean, except for precipitation and stream discharge, which were summed according to calendar year (January-December) and to water year (preceding October-September).

Human impact datasets
To understand how wetland conditions are influenced by human-related impacts in the region, we collected publicly available data on urban development, wetland protection status, and wetland restoration activities for each site (Appendix S2).Watershed development was quantified for each wetland as the percent developed areas in 2016 using the Landsat-derived National Land Cover Database (NLCD; Yang et al., 2018).Developed areas included residential and commercial development with impervious surfaces ranging from 20% to 100%.
We recorded the presence/absence of mouth armoring (jetties, seawalls) and channelization of associated streams and rivers at each site using the coastal armoring dataset (Dare, 2005) and visual inspection of aerial imagery.We used the California Protected Areas Database to determine the level of protection enforced at each wetland site.Public access levels to the protected open spaces are categorized as open access, restricted access, or closed access.We included sites within military land use areas as closed access.
We compiled a dataset of restoration projects through the Southern California Wetlands Recovery Project, the California State Coastal Conservancy, EcoAtlas, and published reports and websites (Appendix S1: Table S2; Appendix S2).For each restoration project, we collected the dates that each project was implemented on the ground, the estimated area to be restored, and the type of restoration activity.We limited the dataset to larger restoration projects that impacted >10% of the wetland site, which pertain mostly to active earth-moving, dredging, and wetland enhancement projects.

Trends in wetland greenness and area
We estimated trends in NDVI and wetland area for each wetland over the study period (1984-2019) using Sen's slope (Sen, 1968).We compared trends with physical site characteristics including latitude, elevation, size, archetype classification, habitat composition, wave climate subregion, and estuary mouth closure.We also compared trends with site characteristics related to human impacts, which included watershed development, coastal armoring, protection status, ownership, and restoration activities, as these are nonclimatic factors that may explain long-term changes in wetland conditions but were not available as continuous  of (A) mean sea level for five regional NOAA Stations, (B) mean monthly significant wave height per wetland site, (C) mean stream discharge for 11 regional USGS gauges, (D) regional means of temperature (black) and precipitation (blue), (E) regional means of drought indices, and (F) regional means of climate indices.Note dual axes for temperature and precipitation data (D) and sea level and wave height reported in the North American Vertical Datum of 1988 (NAVD88).MEI, Multivariate ENSO Index; ONI, Oceanic Niño Index; PDO, Pacific Decadal Oscillation; PDSI, Palmer drought severity; SPEI, Standardized precipitation-evapotranspiration index; SPI, Standardized precipitation index.

Environmental category
time series.For sites that have undergone large restoration projects (Appendix S1: Table S2), we estimated pre-and post-restoration trends.Statistical significance was tested using simple linear regression for quantitative variables, one-way ANOVAs for categorical variables, and paired t test for pre-and post-restoration trends using the base R package "stats" (v1.1-18).

Drivers of wetland greenness and area
We used generalized additive models (GAMs) to describe the complex relationships among annual wetland response (NDVI, wetland area) and environmental drivers.GAMs allow us to model variability between groups (sites, archetypes) in regression relationships, do not require strict assumptions of linearity, and can facilitate several different smooth functions to describe potentially nonlinear relationships, unlike generalized linear models (Wikle et al., 2019;Wood, 2017).We used thin-plate regression splines as the smoothers in our models (Wood, 2003) and penalized models to restrict overfitting of the smooth functions, which can be used to identify linear effects and to select appropriate model covariates (Marra & Wood, 2011).We also applied hierarchical GAMs (HGAMs), an extension of GAMs that help model between-group variability in regression relationships by including a random effect term in the model that accounts for groups within hierarchically structured data (Pedersen et al., 2019).
We fit GAMs to identify common (i.e., regional) drivers to wetland NDVI and area for all sites in the region using the R package "mgcv" (z1.8-36;Wood, 2011Wood, , 2017)).Models were fit to our data using a log-link function under a Gaussian error distribution for the family and basis functions, and restricted maximum likelihood (REML) for the smoothing parameter (Wood, 2011).We compared models based on Akaike information criterion, degrees of freedom, and deviance explained to balance model fit and complexity (Wood et al., 2016).We initially tested all environmental covariates using penalization to aid model variable selection and reduce multicollinearity within driver categories of sea level, waves, temperature, precipitation, stream discharge, drought, and climate (Appendix S1: Figure S3).Ultimately, we selected covariates based on the penalization analysis and our hypotheses on what drives wetland conditions in this region.Once the best model was identified for the region, we then included random effects of sites and archetypes as model terms to test for variability in model intercepts among groups using the HGAM structure (Pedersen et al., 2019).For grouping terms (sites and archetypes), basis functions were set to random effects and basis size (k) was set to match the number of groups.

Long-term changes in wetland greenness and area
Long-term NDVI trends from 1984 to 2019 ranged from −0.0016 to 0.0048 NDVI year −1 for individual wetlands in the region (Figure 2).Over half of the wetlands (18 out of 32) exhibited significant positive trends in NDVI.One wetland, the Santa Ana River Wetlands, exhibited a significant negative trend of −0.0012 NDVI year −1 .The 13 remaining wetlands had nonsignificant changes in NDVI over the 35-year study period.
Wetland area trends indicated significant loss, that is, type conversion of vegetated or mudflat habitats to subtidal, occurring in 15 wetlands.Losses ranged from 0.05 ha year −1 in the Ormond Beach Wetlands (0.3% loss of site area year −1 ) to 3.96 ha year −1 in the Bolsa Chica Lagoon (1.7% loss of site area year −1 ), with an average of 0.80 ha of wetland area being converted to subtidal per year.Significant gains of wetland area occurred in five wetlands, with land area trends ranging from 0.01 ha year −1 (0.1% gain of site area year −1 ) in Las Flores Creek to 0.77 ha year −1 (0.6% gain of site area year −1 ) in San Elijo Lagoon, while 12 wetlands showed no significant changes in wetland area.Overall, wetlands displayed three types of long-term response (Figure 2): greening and gaining wetland (10), greening and losing wetland (16), and browning and losing wetland (6).No wetlands exhibited browning and area gains.
Trends in NDVI were overall not significantly correlated with the site characteristics tested, including latitude, elevation, size, archetype, habitat composition, wave climate subregion, mouth closure, watershed development, protection status, mouth/channel engineering, or restoration activities (Figure 3).Trends in wetland area were significantly correlated with site elevation and size (Figure 4).Mean site elevation was positively correlated with wetland area trends (r 2 = 0.27, p = 0.002), indicating reduced wetland loss in sites with higher mean elevation.The correlation between site size and area trends was negative (r 2 = 0.27, p = 0.002), with larger sites losing more wetland area over 35 years.
Noteworthy differences in the variability of NDVI and wetland area trends among archetypes (Figures 3 and 4) suggest that intermediate estuaries and fragmented river valley estuaries have exhibited a wider range of responses to similar drivers over 35 years.Small creeks and lagoons, however, responded similarly in terms of both NDVI and area trends, with all sites exhibiting greening and gains in area.
Estuarine mouth dynamics and mouth modifications were site characteristics that showed patterns with wetland NDVI and area (Figures 3 and 4).Wave climate subregions reflect varying coastal exposure and setting across the southern California region, and we found that estuaries along steep, terraced gradients and low-exposure coastlines (e.g., Gaviota Coast, SM Bay) had narrower ranges of NDVI and wetland area response.Larger ranges of responses occurred in westerly facing, higher exposure subregions (e.g., Ventura Coast, SP Bay, SD Coast).Sites with predominantly open-mouth closure states showed greater variability in NDVI and area trends and exhibited the majority of browning and wetland loss responses (Figures 3 and 4).Furthermore, sites that were predominantly open due to jettied estuary mouths tended to have more negative trends in NDVI and area (Figures 3 and 4).
Pre-and post-restoration trends for sites with completed projects highlight the contextuality of wetland restoration, as no single pattern emerged in restoration's impact on NDVI or area trends for all restored sites (Figures 3 and 4).The average trend in NDVI across restored sites decreased from pre-restoration (0.000485 ± 0.005 NDVI year −1 ) to post-restoration (−0.000534 ± 0.004 NDVI year −1 ), but differences were not significant (Figure 3; t(9) = 0.31, p = 0.76).Post-restoration NDVI trends increased in six sites, with Sweetwater Marsh exhibiting the most pronounced increase (0.016 NDVI year −1 ) 20 years after restoration, while NDVI trends decreased post-restoration in four sites: Bolsa Chica Lagoon, Santa Ana Wetlands, Carpinteria Salt Marsh, and Upper Newport Bay.Pre-and post-restoration trends in wetland area also did not significantly differ (Figure 4; t(9) = −1.04,p = 0.32), but seven sites exhibited increases in area trends, while three sites exhibited decreases in F I G U R E 2 Regional and site-based trends in wetland conditions in southern California from 1984 to 2019.Regional normalized difference vegetation index (NDVI) trends are categorized as greening, browning, and no or nonsignificant change; Wetland area trends are categorized as loss, gain, and no or nonsignificant change; Wetland response is categorized as greening + gain (green), greening + loss (light green), and browning + loss (brown).Site-based trends show wetland area and NDVI, with nonsignificant NDVI trends indicated by black outlines.Labels correspond to site names indicated in Appendix S1: Table S1, Figure S1.post-restoration: Huntington Beach Wetlands, Carpinteria Salt Marsh, and Upper Newport Bay.The Devereux Lagoon and San Elijo Lagoon restorations were too recent to test for post-restoration trends.Overall, NDVI and area trends were not found to correlate with time since restoration.

Regional environmental conditions and drivers of wetland response
Southern California experienced significant positive trends in mean sea levels and minimum annual air temperatures from 1984 to 2019 (Table 1).Sea levels in the region increased at an average rate of 2.0 mm year −1 (p < 0.0001) and minimum temperatures increased at a rate of 0.025 C year −1 (p < 0.05).No other significant trends in other environmental drivers were found.
The best predictors of regional wetland NDVI and area for all sites were mean sea level, mean wave height, total annual precipitation, annual stream discharge, VPD, and mean annual air temperature, which together explained 47% of the deviance in NDVI but only 27% of deviance in wetland area (Table 2, Model 2).The effective df (edf) for significant covariates indicate near-linear effects of wave height (edf = 0.96) and sea level (edf = 1.24), and more complex effects of stream discharge (edf = 3.54) and VPD (edf = 2.50; Table 2, Model 2, Figure 5).Sea level, precipitation, and mean wave heights showed an increasingly positive effect on NDVI (Figure 5).The mean effect of temperature was nonlinear and nonsignificant but indicated a peak positive effect at 14.5 C, while temperatures above 17 C have a negative effect on NDVI.The effect of stream discharge on NDVI was mostly negative, except for small amounts of inflow (~300 m 3 s −1 ).VPD had a nonlinear relationship with NDVI, with positive effects up to 10.5 hPa and negative effects from 10.5 to 16.5 hPa.
The random effects of archetypes increased the deviance explained in NDVI to 52.7%, but did not substantially change the predicted effects of model covariates with the exception of precipitation (edf = 1.95;Table 2,  Model 3).Grouping at the archetype level therefore indicates similar effects of sea level, wave height, precipitation, stream discharge, and VPD on wetland NDVI within wetland types.However, the random effects of sites had the most significant effect on NDVI, followed by temperature (edf = 0.7), and this model explained 86% of the deviance in NDVI (Table 2, Model 4).
Testing the same covariates in GAMs for wetland area changes showed temperature, wave height, and VPD as T A B L E 2 Comparison of regional models that describe wetland normalized difference vegetation index (NDVI) and area based on top environmental predictors.having the strongest effects (Table 2; Appendix S1: Figure S4).However, the maximum deviance in percent area explained in regional models was only 26.5%.When the random effects of archetypes and sites were included, deviance explained increased to 58% and 95%, respectively, indicating that changes in wetland areas over the period of this study are only partially explained by regional environmental drivers without considering site-specific or contextual impacts at local scales like restoration efforts, nearby infrastructure, or watershed or inlet management.

DISCUSSION
Coastal wetlands in southern California displayed mixed trends in wetland greenness and habitat composition from 1984 to 2019.The majority of wetlands in the region (26 out of 32) showed increased greenness over this period.Greenness increased in wetlands that both lost vegetated marsh and mudflat area (16 wetlands), as well as in wetlands that increased in area (10 wetlands).Although most wetlands showed increased greenness, six wetlands had trends of browning coincident with wetland area loss.
The major drivers of wetland greenness over this study were sea levels, wave height, precipitation, stream discharge, VPD, and air temperature.We expected that the greenness of southern California's wetlands would be impacted by precipitation, similar to other drought-prone Mediterranean ecosystems worldwide (Allen, 2003).In a climate with wet winters, dry summers, and highly variable precipitation, recent low rainfall and high potential evaporation have led to drought conditions in the region (MacDonald, 2007).We found that precipitation and stream discharge showed the largest magnitude of effects on wetland NDVI.As precipitation increases, it has a positive effect on NDVI.However, increasing mean annual stream discharge leads to negative effects on NDVI, which could be due to prolonged inundation, causing marsh die-back or increased sedimentation.Mean annual precipitation increased across the region from the start and end of the study period, yet the overall annual trend showed a decline of −0.16 mm year −1 (Table 1).Stream discharge also showed an overall declining trend of −3.35 m 3 s −1 , and decreases in mean annual discharge over the study period (Table 1).Thus, more precipitation and less stream discharge contribute to positive effects on wetland greenness, highlighting the sensitivity of F I G U R E 5 Regional effects of environmental drivers on annual wetland normalized difference vegetation index for the best explanatory generalized additive model (Table 2, Model 2).VPD, vapor pressure deficit.
semiarid wetlands to factors influencing the timing and magnitude of pulse events affecting wetland hydrology.
Altered precipitation patterns and streamflows can also result in plant stress (Zedler et al., 1986), which may explain the negative effect of increasing stream discharge on NDVI.Changes to precipitation can impact hydrological regimes, soil salinity, and wetland geomorphology over time (Day et al., 2008).In southern California, peak streamflow and precipitation rates are strongly correlated with sediment discharge (Warrick et al., 2015) and urbanization increases peak discharges and runoff volume (Beighley et al., 2003), leading to increased eutrophication on local scales (Howard et al., 2014).In the future, extreme precipitation events in California are expected to intensify, bringing more rain over shorter periods (Dettinger et al., 2011;Huang et al., 2020).The timing and duration of impacts are important to plant response, as shown through pulse disturbance events like fires (Brown et al., 2019) and rapidly elevated sea levels coinciding with El Niño events (Goodman et al., 2018;Harvey et al., 2020) that cause local plant mortality and marsh degradation.
Increasing sea levels and wave heights had increasingly positive effects on NDVI, where NDVI was measured in only vegetated marsh and unvegetated mudflat areas (Figure 5).Sea levels showed significant, positive relationships with NDVI for several sites ranging in size and archetype in the south (Appendix S1: Figure S5), where sea level rise rates are on average 0.93 mm year −1 higher than in the north.The Tijuana River Estuary is one southern site where sea level rise may help combat high sedimentation from river inflow by improving tidal influence (Callaway & Zedler, 2004;Wallace et al., 2005).Increased tidal inundation may also be beneficial by increasing estuarine flushing (Clark & O'Connor, 2019;Largier et al., 2019).Increased flooding and rainfall decrease soil salinity, which can cause increases in productivity and biomass, and lead to changes in marsh species composition over time (Callaway & Sabraw, 1994;Zedler, 1983).
Air temperature and relative humidity play an important role in plant productivity, as shown through the complex nonlinear effects of VPD and temperature on wetland NDVI.Temperature showed a small positive effect on NDVI below ~17 C, but then the effect decreased with higher temperatures.The effect of VPD on NDVI is positive but decreases as VPD nears a threshold of 10.5 hPa, above which the effect on NDVI is negative.VPD increases exponentially with temperature, and plants have limits of tolerance for VPD before negative effects to stomatal conductance and transpiration, which leads to declined plant photosynthesis and growth (Grossiord et al., 2020).VPD has been identified as a major contributor to recent drought-induced mortality in terrestrial plants (Grossiord et al., 2020).In semiarid coastal wetlands, temperature increases can lead to drought conditions and hypersaline soils, which can negatively affect plant community diversity and health (Kelso et al., 2020;Wigginton et al., 2020).
Work from other coastal regions highlights the interacting drivers of temperature and precipitation in wetland conditions.In tidal marshes of the Pacific Northwest, USA, peak biomass was positively correlated with total annual precipitation, while the timing of peak biomass was negatively correlated with average growing season temperature based on 31 years of data (Buffington et al., 2018).From Georgia, USA, wetlands also showed a positive, significant relationship between AGB and total precipitation, minimum temperature, river discharge, and mean sea level over a 28-year period (O'Donnell & Schalles, 2016).In coastal Louisiana, USA, the timing of peak phenology in saline marshes was correlated with air temperature; however, areas of brackish and saline marshes declined from 1984 to 2014 in relation to increasing sea levels, temperatures, and CO 2 (Mo et al., 2019).Findings from these regions emphasize that macroclimate drivers can elicit similar general responses, but that location is important in modeling wetland characteristics (Buffington et al., 2018), that long-term response varies among wetland types (Mo et al., 2019), and that wetland response varies across seasons (O'Donnell & Schalles, 2016).
Wetland area loss was a common trend among 22 of the 32 sites monitored from 1984 to 2019 in this study.Sea levels were identified as the most important driver of wetland area among sites, and overall sea levels in the region have increased by ~2 mm annually.Thus, sea level rise has already been affecting the extent of wetland habitats in specific sites across the region from 1985 to 2019.The Mugu Lagoon and Huntington Beach Wetlands are two examples of non-restored sites that exhibited significant negative trends in wetland area.Annual wetland area at these sites was significantly, negatively correlated with annual sea levels (Appendix S1: Figure S6).Without considering the impacts at individual sites, the combined effect of sea levels, waves, stream discharge, precipitation, and temperature explained 27% of the variability in percent area of wetlands across the entire region.Much of the variability in both wetland area and NDVI was explained by differences among sites.Therefore, regional climate factors play a critical role in the trajectory of wetland change, but site-specific or contextual factors at local scales are needed to further explain the full range of responses we measured in southern California's wetlands.
Local environmental factors, human influences, and physical site characteristics did show patterns with long-term response in wetland area.Site size was positively related to wetland loss, while site elevation had a negative relationship with loss.Our predictions of future sea level rise response across southern California showed a similar pattern of greater expected habitat loss in larger wetlands containing large proportions of vegetated marshes and unvegetated mudflats (Doughty et al., 2018).Large, low-lying wetland archetypes like large lagoons, large river valley estuaries, and intermediate estuaries had some of the worst trends in wetland loss detected in the region as these are more susceptible to sea level increases.Without adequate sediment supply, accretion rates in low-elevation wetlands are unable to keep pace with sea level rise (Thorne et al., 2018).Systems open to the ocean have also shown decreased accretion rates and marsh drowning compared with higher elevation intermediate estuaries (Thorne et al., 2021).
Estuaries with opened inlets, whether natural or engineered, showed greater propensity for negative trends in NDVI and area.Intermediate estuaries, fragmented river valley estuaries, and large lagoons exhibited more variability in greenness and area trends and are some of the most vulnerable archetypes.Intermediate estuaries, or intermittently opening and closing estuaries, are key ecological habitats along California's high-energy coastline (Clark & O'Connor, 2019).Being wave-dominated systems sensitive to fluvial inputs, estuary mouth dynamics, and modification, they are especially geomorphically diverse (McSweeney et al., 2017), which is reflected in the wide range of responses.For example, San Onofre Creek had the highest trend (0.004 NDVI year −1 ) coinciding with detected gains of area, which could indicate reduced tidal and fluvial inputs in a closed system, eventually leading to a transition to nontidal habitats.Conversely, the San Luis Rey Estuary had the lowest trend (−0.002NDVI year −1 ) and relatively consistent habitat composition, suggesting negative effects of ongoing tidal influence to wetland vegetation due to a jettied, predominantly open inlet.Similarly, the Santa Ana Wetlands, a jettied, fragmented river valley estuary, showed a significant negative trend (−0.0012NDVI year −1 ), likely due to this highly modified estuary being disconnected from river discharge and continuously open to tidal changes.Decreased habitat condition and native plant abundance have occurred when natural breaching patterns were altered by bar management (Clark & O'Connor, 2019).Potential resiliency can be broadly characterized by estuary mouth dynamics controlled by coastal exposure and mouth engineering and management.Restoring and adapting intermediate and fragmented estuaries calls for improving hydrological connectivity (SCWRP, 2018), managing environmental flows (Stein et al., 2021), and restoring natural inlet dynamics (Thorne et al., 2021).
Ongoing wetland restoration efforts in the coastal wetlands across southern California play a complex, F I G U R E 6 Observed past trends in wetland loss (this study) compared with predicted future wetland loss (Doughty & Cavanaugh, 2019).Colors signify past wetland response: greening + gain (green), greening + loss (light green), and browning + loss (brown).Labels correspond to site names indicated in Appendix S1: Table S1.
context-dependent role in area trends from 1984 to 2019.Bolsa Chica Lagoon, for example, exhibited the largest trend in wetland loss (0.05 ha year −1 ) detected, but this long-term trend is caused by a rapid conversion of wetland habitats to subtidal over a 3-year period coinciding with a large, earth-moving restoration effort.The Bolsa Chica Wetlands Restoration Project occurred from 2004 to 2007, with the goal of increasing subtidal habitat and restoring full tidal influence to hundreds of hectares (Appendix S1: Table S2).San Elijo Lagoon showed the largest trend in area gains (0.77 ha year −1 ) from 1984 to 2019.In this instance, detected areal increases may indicate long-term degradation to wetland habitats caused by disrupted hydrology and increased sediment loads from surrounding development that led to conversion to nontidal habitats (Sutula et al., 2006).The San Elijo Lagoon Restoration Project began construction in 2017, with the goal of restoring tidal circulation by widening and deepening estuary channels to combat eutrophication and rebalance wetland habitats.Extensive dredging, however, can also have detrimental impacts on vegetated wetland habitats (Kennish, 2001).
Understanding the influence of humans and environmental drivers on coastal wetland health in an urbanized region is complicated by variability and context.By quantifying past trends, we are reminded to prepare for the worst scenario and be aspirational in our goals to mitigate future losses (Figure 6).The important work of protecting and restoring these systems regionwide requires us to account for variability and build a holistic perspective informed by our knowledge of historic, present, and future conditions (Stein et al., 2020).It also requires consistent definitions, goals, and assessment of restoration.Wetland archetypes provide a framework for grouping similar systems to make sense of their diversity and their responses to climate-and human-driven changes (Stein et al., 2020).Regional restoration efforts could benefit from standardized methods to establish baselines of wetland conditions and track future wetland health using remotely sensed data (OPC, 2020; Shuman & Ambrose, 2003).

ACKNOWLEDGMENTS
The views and opinions expressed in this article are solely the authors' and do not reflect the positions, policies, or opinions of the Texas General Land Office or the State of Texas.This work was conducted with funding awarded to Kyle C. Cavanaugh from the National Aeronautics and Space Administration's New Investigator Program (NNX16AN04G).Cheryl L. Doughty was also supported by teaching assistantships through the Department of Geography at University of California, Los Angeles.

F
I G U R E 3 Normalized difference vegetation index (NDVI) trends related to physical site characteristics and human impact variables.Wetland response is categorized as greening + gain (green), greening + loss (light green), and browning + loss (brown).Archetypes include small creeks and lagoons (SCL), intermediate estuaries (IE), large lagoons (LL), large river valley estuaries (LRVE), and fragmented river valley estuaries (FRVE).Wave climate subregions include the Santa Monica (SM) Bay, San Pedro (SP) Bay, and the San Diego (SD) Coast.Boxplots represent data means (midlines), 25th and 75th percentiles (box limits), data values within 1.5 times the interquartile range (whiskers), and outliers.

F
I G U R E 4 Wetland area trends related to physical site characteristics and human impact variables.Wetland response is categorized as greening + gain (green), greening + loss (light green), and browning + loss (brown).Archetypes include small creeks and lagoons (SCL), intermediate estuaries (IE), large lagoons (LL), large river valley estuaries (LRVE), and fragmented river valley estuaries (FRVE).Wave climate subregions include the Santa Monica (SM) Bay, San Pedro (SP) Bay, and the San Diego (SD) Coast.Boxplots represent data means (midlines), 25th and 75th percentiles (box limits), data values within 1.5 times the interquartile range (whiskers), and outliers.