The Altamaha River is one of the largest freshwater sources to the Atlantic Ocean within the United States. It has a significant impact on the Georgia coast and forms a complex system of distributary channels, creeks, and intertidal areas (Di Iorio and Castelao, 2013). The Altamaha estuary experiences semidiurnal tides with an average amplitude of ~2 m; intrusion of saltwater is limited to approximately the lower 20 km of the river (Alber and Sheldon, 1999). The polyhaline portion of the estuary has extensive salt marshes dominated by Spartina alterniflora (Więski and Pennings, 2014).

Our study encompassed S. alterniflora marshes located within the lower Altamaha tidal watershed ( Figure 1 ; light blue), as delineated by the US Geological Survey HUC boundary (USGS HUC 03070106; see https://water.usgs.gov/GIS/huc.html) and an overlay of S. alterniflora communities (green shading) along the central Georgia coast. Recent digital mapping of Georgia’s coastal saltmarsh and brackish marsh at 2 m resolution (Alexander and Hladik, 2015) was resized using ENVI 5.3 (Harris L3Geospatial) to 30 m resolution, in order to match the 30 m pixel resolution of the Landsat 5 imagery used in this study ( Figure 1 ). Within the Altamaha tidal watershed, there are 43,847 pixels of S. alterniflora, 30 m in size, covering ~ 39.46 million m² (3,946 ha).

Spartina alterniflora aboveground biomass estimates were obtained from satellite images as described in O’Donnell and Schalles (2016). Briefly, the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) protocol for atmospheric correction of Landsat 5 and Landsat 8 images was used to estimate Spartina biomass in salt marshes within the study site. The nominal return rate for Landsat imagery is 16 days; however, two overlapping Landsat Path/Rows (16/38 and 17/38) cover the study area (O’Donnell and Schalles, 2016), providing two prospective scenes every 16 days. Many satellite images were discarded because of cloud coverage and tidal inundation, which obscured the canopy reflectance signals of coastal marshes. As described in O’Donnell and Schalles (2016), no images coinciding with a tidal inundation recording above 0.65 m relative to NAVD88 at the Fort Pulaski gauge (NOAA station 8670870, the closest to the study site) (offset by +25 minutes) were used.

The analyses used in this study require equally spaced timeseries. Therefore, from the available satellite imagery, monthly means of Spartina alterniflora aboveground biomass (g/m²) were calculated from the beginning of 1984 to the end of 2018. However, due to cloud coverage, inundation, and lack of Landsat data from late 2011 to early 2013, the availability of suitable images varied widely. Over the 35-year observational period, there were a total of 227 missing monthly biomass estimates out of 420 months (193 existing data points), with data gaps ranging from 1 to 17 months.

To establish a continuous time series of biomass at monthly resolution, we tested several gap-filling approaches including autoregressive integrated moving average (ARIMA) models, kriging, and inpainting. ARIMA models utilize autoregressive and moving average components of time series data to predict (Box et al., 2015) or gap-fill (e.g., Afrifa-Yamoah et al., 2020; Dorich et al., 2020) time series. We used the R package imputeTS (v3.2), which seasonally decomposes time series, fits an ARIMA model to the time series, uses Kalman smoothing to fill the gaps, and then adds the seasonal component again (Moritz and Bartz-Beielstein, 2017). Kriging is a geostatistical approach commonly used for spatial interpolation (Li and Heap, 2008). As suggested by Knotters and Heuvelink (2010) and Lepot et al. (2017), this approach can be adapted to temporal interpolation such as predicting missing values in air temperature time series (Shtiliyanova et al., 2017). In this study, ordinary kriging was performed. Due to the strong seasonality in marsh aboveground biomass, the data were recast in two dimensions, with the time of year (i.e., month) as one axis and the calendar year of an observation as the second axis, similar to the approach of Walter et al. (2013). This 2D kriging method was implemented using the R package gstat (v2.0-9) (Pebesma, 2004). From the available variogram models the exponential model was selected as it produced the best fit for our biomass data. Since satellite data coverage was more complete during winter than summer months due to reduced cloud cover, the seasonal axis was set to start in January, ensuring that the regions of data gaps were well embedded into measured data. For the inpainting approach, we again separated the seasonal pattern from interannual long-term trends. Missing biomass values in the 2D month vs. year coordinates were in-filled by solving a boundary value problem for an elliptic partial differential equation. This was implemented using the MATLAB script inpaint_nans (their spring-metaphor method) by D'Errico (2022).

The performance of gap-filling methods was assessed both quantitatively and qualitatively. The former was measured by the accuracy of predicting existing data points with bootstrapping using 100 Monte Carlo realizations and computing prediction uncertainties by comparing modeled and observed data. The qualitative assessment focused on periods with long data gaps (such as the 17-month period between the last data from Landsat 5 in November 2011 and the first data from Landsat 8 in March 2013) where the lack of data constraints can lead to poor gap-fill performance. Reproducing seasonal patterns was also critical since late Summer/Fall season peaks in biomass are important metrics of salt marsh productivity (Visser et al., 2006; Kirwan et al., 2009).

Monthly records of air temperature, precipitation, sea level and Palmer Drought Severity Index (PDSI: high = wet; low = dry) for 1984-2018 were obtained from the NOAA National Climatic Data Center (NCDC). Monthly mean air temperature and total precipitation data were acquired from the NOAA Sapelo station ( Figure 1 , Lawrimore et al., 2016). PDSI data were retrieved as a monthly time series for the Southeast Climatic Division (NOAA, 2022). In contrast to the biomass data, these datasets have minimal gaps (~10%), which were gap-filled with kriging, as this was the method selected for gap-filling the biomass data (see Results). River discharge data for the Altamaha River were obtained from the USGS Doctortown station because it is the most downstream river gage that covers the time period of our study and provides an accurate measure of all of the water that enters the estuary ( Figure 1 , U.S. Geological Survey, 2016). Nutrient concentrations (NO_x and PO₄ ^3–) were obtained from the USGS Altamaha River Near Gardi station which is approximately 30 miles upstream from the Altamaha Estuary ( Figure 1 , U.S. Geological Survey, 2016b). We used nutrient concentrations as a proxy for the composition of the water flooding the marsh to which the plants are exposed. The small number of missing values (<10%) in the nutrient dataset were also gap-filled by kriging. Monthly mean sea level data and tidal projections were obtained from the NOAA Fort Pulaski station ( Figure 1 ; NOAA 2022a; NOAA 2022b).

The continuous wavelet transform (CWT) was used to analyze the time series in time/frequency space. We selected Morlet wavelets for their balance of time and frequency localization and facilitation of detection of time-dependent amplitude and phase (Lau and Weng, 1995; Grinsted et al., 2004). To investigate local correlations between CWTs of environmental variables and biomass, coherence patterns in time/frequency space were computed and phase shifts were quantified. Statistical significance of the results was estimated using Monte Carlo randomizations. Both wavelet and wavelet coherence analyses were implemented using the biwavelet (v0.20.21) R package (Grinsted et al., 2004).

Ensemble empirical mode decomposition (EEMD) was used to decompose the time series into intrinsic mode functions (IMF) (Wu & Huang, 2009). This data-driven method analyzes nonlinear and nonstationary processes by breaking down complex time series into components that are based on the observed minima and maxima in the signal. By sequentially subtracting these modes from the original signal, the time series is effectively decomposed into signals of increasingly lower frequencies (modes). To address the issue of the occurrence of similar frequency signals in multiple modes (mode mixing), IMFs are computed as the ensemble average of multiple stochastic realizations of empirical mode decomposition (Huang et al., 1998), after addition of white noise to the original data. Noise was set to 20% of the standard deviation as suggested by Wu and Huang (2009), and the MATLAB implementation of Yang et al. (2018) was used with Z-score normalization.

To explore univariate causal connections between environmental variables and salt marsh biomass, we used convergent cross mapping (CCM; Sugihara et al., 2012). This empirical dynamic modeling approach builds on Takens’ theorem (Takens, 1981), creating state-space reconstructions from lagged time series. The dimensionality of the state space and prediction time lags between biomass and other variables were determined as the ones that maximize prediction skill (Ye et al., 2015). A causal connection was then inferred from the information embedded in the response variable (i.e., biomass) that leads to higher prediction skills with an increasing length of the time series (library size used in the reconstruction of the manifold). Prediction skills produced by causal variables were expected to be higher than regular cross correlation between the two variables (Bonotto et al., 2022). This method was implemented using the R package rEDM (v1.9.3).

To avoid potential false causal implication resulting from seasonality in two variables, we compared our analysis of the observed data with outcomes using randomized seasonal time series (Deyle et al., 2016). One thousand surrogate time series were generated with a seasonal pattern based on the multi-year monthly averages and adding randomly shuffled residuals to each point. These surrogates were tested against biomass in CCM, and the resulting prediction skills were obtained. By comparing those with the prediction skill computed for the real data, the confounding effect of seasonality was quantified.

A strong seasonal pattern was apparent in our biomass data ( Figure S1 , top panel). Notably, imagery (and derived biomass estimates) was most limited during summer months because of frequent cloud coverage during months with high humidity, precipitation, and temperature (bottom panel in Figure S1 ; see also O’Donnell and Schalles, 2016), which led to larger uncertainties in the biomass estimates during peak season. Additionally, interannual differences in biomass led to larger variability during periods of substantial growth or decay (spring green-up and fall senescence), contributing to larger uncertainties during those parts of the year.

The different gap-filling methods for the monthly biomass data showed similar performance. ARIMA/Kalman performed best when data gaps were small, but yielded the lowest accuracy when the number of missing data points was high (i.e. > 50%) ( Figure S2 ). The two non-traditional temporal interpolation methods, inpainting and kriging, with years and time of year as major axes, performed comparably to or better than the ARIMA/Kalman approach, depending on the amount of missing data. Both yielded similar results and produced realistic seasonal peaks. Given the similar gap-filling performance, datasets filled by kriging were used in this study because of its statistical basis and more extensive literature support than inpainting (Shtiliyanova et al., 2017).

Monthly averaged biomass ranged from 254 g/m² to 2650 g/m², and gap-filling could recover missing peaks in biomass data ( Figure 2 ). The largest gap between the end of 2011 and Spring 2013 was caused by the interval between the end of Landsat 5 and the start of Landsat 8 operations.

We used both wavelet analysis and EEMD to identify periodicities in the data. As described below, these yielded similar results.

Wavelet analysis of biomass identified significant periodicity at the episodic (0-4 mo), multi-month (4-8 mo), annual (8-16 mo), and multi-year (64-128 mo) scales ( Figure 3A ). At the sub-annual scale (0-8 mo), intermittent signals of short duration were observed between the late 1980s and late 1990s. There was a strong annual band between 1984 and 2009, after which the annual periodicity became weaker for the remainder of the time series ( Figure 3A ). Multi-year signals were apparent for the period of 1995-2006 (the only years not influenced by edge effects at these long periods). Similar patterns were observed using EEMD, with sub-annual cycles captured in IMF 1, the dominant annual cycle captured in IMFs 2 and 3, and the multi-year signal in IMF 5, with the rest of the IMFs contributing less to the signal ( Figure 4 , left). The EEMD residuals (original signal with all the modes subtracted) showed an overall increase of biomass over time, in particular prior to 2000 ( Figure 4 , top right).

Monthly precipitation on Sapelo Island showed annual patterns, but the signal was discontinuous throughout the time series ( Figure 3B ). It also showed multiple intermittent episodic (i.e., <4 months) signals. EEMD of precipitation also showed these high frequency signals, with the largest signal contributions by intrinsic modes 1 and 2 ( Figure 4 , left). In line with the absence of low frequency signals in the wavelet analysis, intrinsic modes representing variations over longer timescales made a small contribution to the observed precipitation signal ( Figures 4 , S5 left).

Air temperature showed the expected strong and consistent annual pattern ( Figure 3C ). Additionally, signals with longer periodicities (i.e., >1 year) were weaker than the average signal intensities observed in other variables. This pure seasonality was also evident in the EEMD with a dominant IMF2 that exhibited a minimal difference between maximum and average range, and the remaining IMFs contributing considerably less to the observed temperature signal ( Figures 4 , S4 right). The EEMD residual showed an increasing trend of approximately 2°C over 35 years ( Figure 4 , right).

River discharge showed time-frequency patterns similar to precipitation, but with a stronger annual signal and some discontinuous sub-annual signals ( Figure 3D ); this pattern was also apparent in the larger amplitudes in IMFs 1, 2 and 3 ( Figure 4 , left). Discharge also showed higher power at lower frequencies ( Figure 3D ) compared to precipitation.

Annual and sub-annual periodicities in sea level were identified by both CWT and EEMD ( Figures 3E , 4 ); however, these periodicities were intermittent, and most pronounced in 2000 - 2013. The residual in the EEMD analysis revealed the rising trend in sea level during the study period of approximately 0.15 m over the last 3 decades ( Figure 4 , right). Additionally, CWT for NOAA tidal projections showed high seasonality ( Figure S3 ).

The drought index showed periodicities on the order of 8-20 months in the wavelet analysis from the 1980s until around 2000 ( Figure 3F ). In the second half of the time series a shift to a longer period (~5 yr) was observed. EEMD also demonstrated the slower variations of the drought index. It revealed a shift to longer period signals in the second half of the data set with decreased amplitude in IMF2 and increased contributions of IMFs 4 and 5 ( Figure 3 , S6 right). Additionally, the long-term trend showed intensifying drought conditions until the mid-2000s, which aligned with the decrease in precipitation and continuous increase in temperature as underlying variables ( Figure 4 , right). This alignment between environmental variables might have a synergistic effect on biomass; however, since our methods do not consider combined impact of variables on biomass, we could not quantify such effects.

Both riverine NO_x and PO₄ ^3– concentrations showed intermittent seasonality and sub-annual signals ( Figures 3G, H ), which resulted in larger differences between maximum and average ranges of modes 1 and 2 ( Figures 4 , S7 ).

To compare the patterns in hydrological, climatic and nutrient variables with those in the biomass, we analyzed the coherence of their respective wavelets.

The strongest coherence between biomass and precipitation occurred in the annual band, reflecting the seasonality of both biomass and precipitation. Coherence was not persistent over time, and the two variables were mostly in phase ( Figure 5A ).

Temperature and biomass showed coherence in the annual band, which was significant for the period 1984-2014 and weaker afterwards ( Figure 5B ). Phase differences were around 0-3 months, with temperature leading biomass. This in-phase relationship reflected the strong seasonality in both signals, as increases and decreases were aligned. Additionally, longer term signals showed coherence around the 8-year period with a 6.5-year lag between temperature and biomass which was partially affected by the cone of influence.

The wavelet coherence of river discharge with biomass was high in the annual band, with most of the significant values between 1984-2016 ( Figure 5C ). The approximately 6-month lag between discharge and biomass is consistent with the observation that river discharge peaks in early spring and biomass peaks in early fall. After 2008, multi-year patterns (i.e., 3-6 years) of river discharge correlated with biomass, with discharge leading biomass by several months. Note though that the multi-year coherence is affected by the cone of influence and hence has limited support.

Coherence of sea level with biomass showed roughly in-phase annual patterns ( Figure 5D ). From 1984-2001, sea level showed coherence with biomass in longer period signals (~3 years). In these regions, sea level led biomass by approximately 2 years.

Drought index and biomass did not show an annual band in wavelet coherence ( Figure 5E ). The only significant patterns occurred during a 5-year period around 2009 in which PDSI led biomass by less than a year, suggesting that a wet (dry) fall period was followed by high (low) biomass the next year. The drought index was also strongly coherent and in-phase with river discharge over longer periodicities (not shown).

Nutrients generally showed intermittent coherence with biomass. NO_x only had two interpretable short significant zones (annual and 5-year periods), only one of which was outside the cone of influence. This was around the mid-1990s in the annual band, where NO_x led biomass by ~1.5 months. PO₄ ^3– had three similar zones in the annual band, with a slight out-of-phase relationship with biomass. After the mid-2000s, PO₄ ^3– showed out-of-phase coherence with biomass around 3- and 5-year periods.

To assess if the observed coherences between environmental factors and biomass are indicative of a causal connection, the time series observations were analyzed using convergent cross mapping (CCM). Results showed that for air temperature, river discharge, sea level, drought index and river nutrient concentrations, the prediction skill both increased with increasing library size and exceeded their linear cross correlation ( Figure 6 ). Thus, these variables were considered causal to biomass. CCM for precipitation did not yield better prediction skill than maximum lagged cross correlation (with up to 12 months of lag) and the use of longer time series did not steadily improve the predictive power. This indicated the lack of a causal connection.

Tests with surrogate signals were performed to investigate if the apparent causal connections were a result of the inherent seasonal cycles in most of the variables. The results indicated that temperature, river discharge, and sea level data performed significantly better than their seasonal surrogates (p-value <0.05). Similarly, and unsurprising given their weak annual bands ( Figures 5E–G ), the causal connections identified from the drought index and nutrients (NO_x and PO₄ ^3–) to biomass were not due to seasonality. CCM for precipitation did not exhibit prediction skills significantly better than its surrogate, which further supported the lack of causal connection between precipitation and biomass.

To put these observations into context, the patterns in environmental factors potentially affecting plant growth are identified and discussed, and then analyzed regarding their correlation in time-frequency space with that of aboveground biomass using wavelet coherence.

Short-term precipitation signals observed in CWT and captured in the first intrinsic mode reflect the episodic nature of precipitation and fronts generated by storm events (Srock and Bosart, 2009). The wavelet analysis did not reveal significant patterns at longer timescales. However, previous studies have shown that southeastern US precipitation depends on large-scale climate patterns such as the North Atlantic Oscillation, Bermuda High Index, El Niño/Southern Oscillation, and the Pacific Decadal Oscillation (Sheldon and Burd, 2014). Among these, the Bermuda High Index is most influential, especially during drought events when it is negative, pointing to westward migrating high pressure areas that block storm impact on the east coast (Sheldon and Alber, 2013). However, important large-scale oscillations often vary intermittently over longer periods and variable periodicities (e.g., the North Atlantic Oscillation; Markovic and Koch, 2005) and hence may not appear as a consistent pattern in the wavelet scalograms ( Figure 3B ).

Our analyses indicate an alignment of annual peaks in precipitation and biomass ( Figure 5A ). This might reflect that precipitation can be beneficial for biomass as it freshens the marsh environment. Additionally, analysis of precipitation and CO₂ fluxes in a salt marsh (Chu et al., 2021) showed high wavelet coherence in the annual band and revealed the importance of early growing season precipitation. However, precipitation can also act as a stressor. For example, Hanson et al. (2016) have demonstrated the effect of precipitation patterns by showing that S. alterniflora biomass under ambient daily precipitation exceeded that subjected to biweekly storms (with similar amount of total rain) or drought conditions. This indicates that the nature of precipitation events is more important than total amount of precipitation and therefore extreme events might weaken the seasonal precipitation/biomass coherence.

The temperature signal was strongly and consistently seasonal ( Figure 3C ) and EEMD showed an increase of approximately 2°C over 35 years ( Figure 4 , right). This is consistent with long-term warming in the state of Georgia of approximately 3°F since the late 1970s and an overall trend of the Georgia coast warming faster than the state as a whole (Frankson et al., 2022; U.S. EPA, 2022).

Biomass peaks aligned with the temperature peaks, reflecting both of their seasonal cycles. Both biomass and temperature also exhibited an increase over the study period, which suggests a positive relationship between temperature and plant growth. However, the long-term trend in the biomass data indicated a maximum around 2008 while temperatures continue to rise ( Figure 4 ). This points to a potential decoupling between these two variables, which can also be seen in the lower power of the wavelet coherence in the annual band after 2015 ( Figure 5B ). A positive relationship between temperature and plant aboveground biomass may not persist if temperatures increase to levels that decrease plant productivity (35°C, Giurgevich and Dunn, 1979). Such negative impacts are more likely to occur in more productive southern salt marshes studied here, which are closer to their optimum temperatures than northern Atlantic coast marshes (Kirwan et al., 2009).

Discharge showed both annual and longer-term patterns. The annual patterns reflected the historical spring maximum (Alber and Sheldon, 1999). Over longer time scales, there was a drop in discharge between 2000 and 2012 ( Figure 4 , right). This may reflect the potential combined impact of climatic and anthropogenic changes (Takagi et al., 2017). We observed higher power in lower frequency signals for river discharge compared to precipitation ( Figures 3B–D ). This is likely due to the fact that river discharge captures and integrates local effects across the entire watershed, whereas our precipitation data was from a single station on the coast of Georgia. Furthermore, seasonality in evapotranspiration rates may decouple precipitation and discharge on sub-annual timescales.

Coherence of biomass and river discharge shows that with some lag, an increase in discharge leads to an increase in biomass. This can be explained by the effect of discharge on estuarine mixing and the salinity and chemical composition of water that floods the intertidal marshes. The salinity of this water is mainly controlled by the Altamaha River which causes freshening of this estuarine system when discharge is increased (Di Iorio and Castelao, 2013). Consequently, it also affects the porewater salinity, which is an important driver of S. alterniflora production (Odum, 1988; Więski and Pennings, 2014; Miklesh and Meile, 2018). This aligns well with the observed stronger coherence of biomass with discharge than with precipitation.

Sea level showed annual and sub-annual signals as well as longer-term trends. Much of the annual and sub-annual signals in sea level can be attributed to tides (e.g., perigean spring tides), which can be seen by comparing the analyses of sea level with/without removal of the projected tides ( Figures S3 , 3E ). However, discontinuities in significance levels in these bands indicate cases where discrete weather events overwhelmed tidal effects. For example, Andres et al. (2013) have shown that local forcings such as along-shelf wind stress caused interannual variability in sea level along the US east coast. These types of local forcings are highly variable; therefore, they manifest themselves as intermittent pulses and disturbances in the time series, but do not generate persistent signals in the scalograms. The long-term increase in sea level indicated by the EEMD residual was 0.15 m over the last 3 decades ( Figure 4 ; i.e., 5 mm/year) which was comparable to the increase of 3.3 mm/year from 1983 to 2001 reported by Sweet et al. (2022) and showed acceleration in the last 2 decades.

Sea level had considerable seasonal, close to in-phase coherence with biomass, indicating that biomass increases with increased sea level. This positive impact of sea level was reported in previous studies (Więski and Pennings, 2014; O’Donnell and Schalles, 2016); however, this effect was shown to be sensitive to depth and duration of the flooding and can reverse in the long run with certain thresholds exceeded (Morris et al., 2013). This implication was also similar to the findings of Wei et al. (2020), who have shown lagged seasonal and multi-day coherence of tide height and CO₂ fluxes in a salt marsh, and Souza et al. (2022) who documented both in and out of phase coherence between water level and CO₂ fluxes on semi-diurnal (in phase) and multi-day (out of phase) timescales. Our analysis also revealed long period coherences of biomass with temperature (~ 8 years) and sea level (~3-4 years), with lags on the order of >6 years for temperature and approximately 2 years for sea level ( Figure 5 ). These lags are long compared to the Spartina growth cycle and therefore may not reflect biological interactions with biomass.

Periodicities in the drought index shifted from shorter (8-20 months) to longer (5 yr) around 2000. This may reflect recurring extended drought periods, which decreased the impact of seasonality in the PDSI signal ( Figure 3 ; Figure 4 , IMFs 4,5). Long-term patterns in drought and river discharge were similar, highlighting their regional and integrative nature ( Figure 3 ). The correlation between biomass and PDSI as reflected in their wavelet coherence points to the potential impact of droughts. Droughts are known to be associated with large-scale dieback (Alber et al., 2008), and there have been several such events on the GA coast since 2000. Such effects may not be immediate but propagate through feedbacks between above- and below-ground biomass which are tightly connected (Schubauer and Hopkinson, 1984; O’Connell et al., 2021). Note that size-class partitioning of central Georgia Spartina alterniflora (O’Donnell and Schalles, 2016; Zinnert et al., 2021) showed that high-marsh, short form Spartina (canopy < 50 cm) and mid-marsh medium form (canopy 50 – 100 cm) had less resistance to drought-induced biomass declines than low-marsh tall form (canopy > 100 cm) during the 5-yr severe drought event centered around 2000. The shorter forms of Spartina live in naturally stressed, higher pore water salinity areas of the elevation-graded marsh platforms and have lower average above- and below-ground biomass.

Nutrient concentrations in the Altamaha River are affected by increasing anthropogenic inputs resulting from agriculture, livestock, and changes in population density modulated by in-stream processes (i.e., biological/benthic uptake/release, ad-/de-sorption, precipitation/dissolution, dilution) (Schaefer and Alber, 2007; Takagi et al., 2017). Nutrient concentrations were not highly sensitive to discharge (NO₃ ^–: non-significant, DIP: slightly negative significant relationship between discharge and concentration) (Takagi et al., 2017). This could explain the weaker annual band in nutrients compared to discharge ( Figure 3 ). Under drought conditions, short-term nutrient pulses can occur when storms/heavy rains cause flushing of retained nutrients from the terrestrial environment before it is diluted by a high flow event afterwards ( Figure 4 , left, IMFs 1 and 2; Whitehead et al., 2009).

Nutrients showed intermittent coherence patterns with biomass mostly similar to coherence patterns of PDSI and discharge. However, differences in phase lags point to differences in underlying processes ( Figure 5 ).

The causal relationships identified in this study generally align well with observed coherences between environmental variables and saltmarsh biomass described above. Analyses showed that temperature, sea level, discharge, drought, and nutrients were causal to biomass while local precipitation was not.

Regression analyses from previous studies are consistent with our findings. Więski and Pennings (2014), working with annual data between 2000 and 2011, showed that river discharge and sea level best predict salt marsh productivity in salt marshes around the Altamaha River estuary. This is likely due to their impacts on porewater salinity and plant stress, which can be the potential mechanism for the causal connection that we observed in CCM for river discharge and sea level. In this analysis precipitation had only limited correlation with production. O’Donnell and Schalles (2016) came to similar conclusions based on an analysis of Landsat 5 data between 1984 and 2011. They found that river discharge, total precipitation, minimum temperature, and mean sea level were important predictors of aboveground Spartina alterniflora on the Central Georgia Coast. This generally agrees with our finding except for precipitation. Although we observed a temporary seasonal coherence with biomass, a causal relationship was not supported by our analysis. Additionally, PDSI was reported as a predictor in combination with other variables only during limited time periods; our results suggest that drought index is indeed causal to biomass. It has been suggested that the effect of drought on biomass can be through biogeochemical factors (salinity, pH, Alber et al., 2008) and can also intensify ecological controls (e.g., grazing pressure, Silliman et al., 2005). We note that the coherence between the biomass and drought is only strong during the second half of the observational period when drought conditions lasted longer, potentially suggesting that there is a minimum duration required to impact the marsh vegetation.

Spartina alterniflora has been shown to be sensitive to N supply (Morris et al., 2013). Although N is the main limiting nutrient in salt marsh ecosystems, P can also be limiting, as low PO₄ ^3– availability can restrict nutrient replenishment by microbial organisms, in turn impacting carbon fixation in marsh plants (Sundareshwar et al., 2003; Rolando et al., 2022). Supporting these findings, our analysis showed that both nutrients are causal to salt marsh biomass. However, nutrient concentrations on the marsh will differ from those measured at the monitoring station in the freshwater reaches of the Altamaha due to mixing with oceanic water, in-stream processing of nutrients, and local sources and sinks such as those due to nitrogen fixation, nitrogen uptake and denitrification. Considering dilution of riverine nutrients resulted in different causal strength values but still shows the causal connections of nutrients to the marsh biomass ( Figure S10 ). Additionally, the causality of the nutrients was sensitive to the gap-filling procedure for the biomass data ( Figure S11 ). These findings emphasize the need for more detailed experimental studies of the causal relationship between nutrients and biomass beyond the analysis of time series.