Topological and Morphological Controls on Morphodynamics of Salt Marsh Interiors

: Salt marshes are important coastal environments and provide multiple beneﬁts to society. They are considered to be declining in extent globally, including on the UK east coast. The dynamics and characteristics of interior parts of salt marsh systems are spatially variable and can fundamentally affect biotic distributions and the way in which the landscape delivers ecosystem services. It is therefore important to understand, and be able to predict, how these landscape conﬁgurations may evolve over time and where the greatest dynamism will occur. This study estimates morphodynamic changes in salt marsh areas for a regional domain over a multi-decadal timescale. We demonstrate at a landscape scale that relationships exist between the topology and morphology of a salt marsh and changes in its condition over time. We present an inherently scalable satellite-derived measure of change in marsh platform integrity that allows the monitoring of changes in marsh condition. We then demonstrate that easily derived geospatial and morphometric parameters can be used to determine the probability of marsh degradation. We draw comparisons with previous work conducted on the east coast of the USA, ﬁnding differences in marsh responses according to their position within the wider coastal system between the two regions, but relatively consistent in relation to the within-marsh situation. We describe the sub-pixel-scale marsh morphometry using a morphological segmentation algorithm applied to 25 cm-resolution maps of vegetated marsh surface. We also ﬁnd strong relationships between morphometric indices and change in marsh platform integrity which allow for the inference of past dynamism but also suggest that current morphology may be predictive of future change. We thus provide insight into the factors governing marsh degradation that will assist the anticipation of adverse changes to the attributes and functions of these critical coastal environments and inform ongoing ecogeomorphic modelling developments.


Introduction
Salt marshes represent a major component of low-lying sedimentary coastal systems and occur across the world [1]. Over recent decades, salt marshes have attracted increasing attention, with much research being focused on the services and functions they provide [2]. Salt marshes exhibit an extremely high biodiversity [3] and primary productivity [4]. They attenuate wave energy and contribute significantly to the protection provided by natural foreshores from high waves and water levels threatening coastal communities [5,6]. Marshes are a sink for atmospheric carbon [7][8][9], while providing a habitat for many endangered or threatened species and nursery areas for commercial fish stock species [10,11]. They also have cultural value as areas for recreation and tourism. In the UK, over the last 150 years map and aerial imagery suggest an expansion of marsh areas in northern England, as against areal loss in the south, attributed to regional variations in sediment supply [12]. It has, however, proven difficult to ascertain patterns of marsh areal change and controlling factors over the latter half of the twentieth century. The East Anglian coast ( Figure 1) is thought to be a region with high rates of wetland loss (e.g., [13]), but in reality rates and types of marsh loss have exhibited great spatial variability (e.g., [14]). Furthermore, achieving precise estimates of changes in marsh area has been shown to be challenging by the few large-scale UK inventories attempted [15][16][17].
The position of a marsh system within a broader context of the coastal zone exerts controls on factors such as sediment or nutrient import and export, tidal flushing and residence times, and forces exerted by tidal currents. The connectivity between intertidal wetland areas and offshore deep channel zones is crucial to water, sediment, and nutrient exchange and thus to the morphological evolution of marshes [18].
The position of a marsh parcel within the surrounding system, and the position of a point on the marsh within individual parcels, modulate local-scale responses. There has long been recognition that elevation-sedimentation relationships vary with scale. For example, single-tide sediment deposition decreases with distance from tidal channels, while surface elevation provides an important marsh-wide control over annual to decadal timescales [19]. Kearney and Rogers [20] previously used logistic regression to predict internal platform integrity changes in marshes at a regional scale in Chesapeake and Delaware Bays, USA. They demonstrated empirical relationships between changes in marsh surface condition and factors such as distance up-estuary and position within a marsh parcel. The marsh systems on the UK East coast have a different context (tidal range, position within tidal frame, sedimentology, vegetation community) to those studied by Kearney and Rogers [20] but we are able to draw comparisons between the topological relationships presented here and their findings.
The morphology of the marsh, described by the spatial distribution of landscape units (at a scale of metres) such as vegetated platforms, salt pans, creeks, and large channels, is the integrative result of historic morphodynamics [21]. Morphology is also thought to exert a control on future changes through biogeomorphic feedback [22], while the interactions between topography and hydrodynamic forces have been extensively explored from a numerical perspective [23]. From an empirical perspective, a relationship between the functional form of marsh margins and erosion rates has been demonstrated [24].
Of particular, and most immediate, interest for landscape management are the loci of greatest dynamism; the most important to understand when considering future ecosystem service provision are those exhibiting (or likely to exhibit) erosion. This study aims to assess the decadal morphodynamics of salt marsh systems on the east coast of the UK, evaluate the role of topological and morphological factors in determining the observed changes, and provide understanding of the controls on salt marsh morphological evolution to support ecosystem management and the development of models that combine ecological and physical functioning. We present landscape-scale statistical models relating such changes to easily derived spatial parameters at scales from the regional (tens of kilometres) to local (metres). Such understanding will help with the prediction of locations likely to exhibit degradation accompanied by concomitant losses of ecosystem services. This will thereby facilitate targeted management interventions to protect the ecological and physical functioning of these important coastal ecosystems.

Study Area
The coastline of East Anglia, UK, is bounded by the Humber estuary to the north and the Thames estuary to the south ( Figure 1). The region is, in many parts, densely populated. The population in the coastal districts of Suffolk and Norfolk exceeds 600,000 (Office for National Statistics GB 2016 census data-www.ons.gov.uk/census accessed on 23 October 2018). A considerable amount of variability occurs along this stretch of coast in terms of the hydrodynamic and sedimentary contexts. The mean spring tidal range (MSTR) varies from 6.18 m at Immingham on the Humber to a minimum of 1.94 m at Lowestoft, Suffolk, before increasing again further south (ntslf.org accessed on 10 December 2020).  [17] and indicating locations referenced in the text. Blue dashed lines denote the boundaries of ISCE units after [25]. Red dots denote standard ports with tide gauge records.
Pethick and Leggett [25] partitioned the region into three Integrated Scale Coastal Evolution (ISCE) units. The Northern ISCE includes the eroding glacial cliffs north of the Humber, which provide sediment inputs for the sandy shorelines in Lincolnshire and the infilling embayment of The Wash. The ISCE includes the spit and barrier island system of the north Norfolk coast as far east as Cromer. The second ISCE is dominated by cliffs composed of glacial sands and gravels and lies between Cromer and Thorpeness. The third, southerly, unit comprises numerous estuaries and inlets sitting within a large embayment and characterised by silt or clay sediments [26]. The region contains some 14,406 ha of salt marsh [17] of diverse character and setting; open coast, embayment, back-barrier, and estuarine marshes are represented [27]. Suffolk and Essex have both experienced a net loss of marsh area in the second half of the twentieth century [15] while The Wash embayment, between Norfolk and Lincolnshire, represents a long-term sediment sink [28] and continues to infill [26], with attendant increases in marsh area [24]. At the regional scale, none of these systems can be thought of as 'sediment starved'.

Morphological Change
We address controls on morphodynamic change within marsh interiors comprising the complex morphologies of vegetated surfaces, creeks, pools, and pans that lie landward of the seaward margin of the marsh as a whole. We use satellite imagery to estimate changes in vegetation distributions within these areas, from which we infer morphological changes. This inference is based on well-established elevational controls on intertidal vegetation establishment [29], which have been thoroughly ground-referenced for NW Europe by Suchrow and Jensen [30]. In inferring morphological change, we assume that all surfaces of sufficient elevation to support vegetation become colonised, while those that are too low do not.
Changes in the extent of vegetation within marsh interiors were estimated using a modified trend analysis of the Landsat archive, which represents the longest appropriate satellite time series available (1984-present). Imagery dating between 1985 and 2016 was used in this study. The exact time period over which the metric is calculated will vary slightly between pixels due to the different acquisition dates of imagery over certain areas and cloud cover precluding the use of some pixels (see Supplementary Material (S1)).
A metric was derived, denoted δP veg , that can be interpreted as representing the percentage change in vegetation cover within any given pixel over the timeframe of the satellite observations. As such, it reflects the temporal variations in the areal unvegetated-vegetated marsh ratio (UVVR) within each 30 m by 30 m pixel. The UVVR could be considered a geomorphic metric that has itself been shown to be related to marsh vulnerability [31]. The methodological workflow we use is summarised in Figure 2 and is detailed extensively in Appendix A. In Google Earth Engine [32], the Normalised Difference Vegetation Index (NDVI), a proxy for the amount of chlorophyll (and therefore vegetation) present, was computed for all scenes within the Landsat archive for the study region up to 2016. Conceptually, the NDVI for a pixel in a satellite image can be considered to be a function of the percentage of that pixel covered by vegetation and the nature of the vegetation within that pixel. A change in the NDVI reflects a change in the percentage vegetation cover within the pixel, with the signal potentially being modulated by any attendant changes in community composition (and therefore spectral signature) and the vigour of the vegetation present. Lopes et al. [33] evaluted a number of vegetation indices for the monitoring of salt marsh extent and condition in Portugal based on the Landsat archive and concluded that NDVI performed best, with a seasonally varying goodness-of-fit [34] typically exceeding 0.9 at one location and 0.75 at another, where a perfect classification would result in a value of 1. NDVI values were cross-calibrated to account for differences in radiomemtric response between different Landsat sensors before a linear trend was computed over the time series for each pixel. A relationship between the coefficients of the trend and the mean NDVI over the time series was found, which was removed by taking the residuals of the regression fit between the slope of the trend line and the mean NDVI (R 2 = 0.22, p ≤ 0.001). The resulting residuals were standardised and calibrated against 25cm Environment Agency aerial photography from 1992 and 2013/14. Calibration was achieved by the visual assessment of the aerial photography subsets extracted for a random sample of satellite pixels (n = 83) stratified to represent the full range of observed standardised residual trend coefficents. For each satellite pixel area, 50 points were randomly distributed and each was manually classified as either coinciding with vegetated or unvegetated surfaces in both the 1992 and 2013/14 aerial photography. For each year, the percentage of satellite pixel area that was vegetated in the aerial photography was then calculated based on the counts in each surface cover class. A strong relationship was found between the satellite-derived residuals and the photography-based estimates of change in vegetation cover (R 2 = 0.87, p ≤ 0.001), implying that changes in vegetation extent dominate changes in δP veg . The calibration sample was drawn from the entire study region. δP veg varies in the range −100% to +100%, where −100% represents a change from fully vegetated to fully bare and +100% the opposite. The method was validated against a sample of pixels representative of the full range of calibrated estimates (n = 100) but drawn from a local subset area (Hamford Water, Essex) to obviate the effect of any large-scale (e.g., latitudinal) signal that may have been present within the calibration. The δP veg metric performed well, with an RMSE of 11.9% of full scale (±100%).

Topological and Morphological Metrics
Two metrics were derived to describe the topology of a location within a marsh and two morphological metrics were developed to describe the distribution of landscape units within a single Landsat pixel area.

Geomorphic Setting
The geomorphic setting of the marsh, meaning its context within the wider coastal system, was represented as a cost function, denoted Cost Distance (CD), describing how difficult it would be for a parcel of water originating offshore to reach any given location within a marsh. Thus, for example, interior marsh areas towards the head of estuaries are more 'costly' to reach (less well connected to the offshore) than the seaward margins of open coast marshes. Without a full simulation of tidal exchanges against which to calibrate flow pathways, the CD metric is not expected to behave in an isomorphic manner with actual tidal exchanges. Rather, CD values are internally consistent and represent a scale of relative connectivity within the study domain.
The area between the 10 km offshore limit and the land was represented on a 10 m grid and divided into four zones. For ease, we denote these 'subtidal', 'intertidal', 'supratidal', and 'terrestrial', while recognising that they do not conform strictly to these descriptions. Intertidal, for example, would ordinarily refer to elevations between the lowest and highest tides experienced. In the UK, this range would extend well below 0 m Ordnance Datum Newlyn (ODN), which approximates mean sea level but is used here as the lower bound of the 'intertidal' zone (see below). The landward limit was defined by the UK Environment Agency Second Generation Shoreline Management Plan as segments vector layer(SMP2), which typically reflects the line of engineered defences. Costs to traverse each cell of the grid were based on distances, with the cells in each zone given different weighting factors. Areas below 0 m ODN were denoted 'subtidal' and were assigned a cost according to their euclidean distance in metres from the offshore boundary.
An 'intertidal' zone was defined as areas between 0 m ODN and the level of the highest astronomical tide (HAT). Elevations were derived from the Environment Agency 2 m-resolution LiDAR composite product (2008) resampled to the 10 m grid. Since spatially resolved data describing HAT were not available, modelled MSTR values [35] were extrapolated to intersect the shoreline and an 0.7× MSTR was used to approximate the level of HAT. Our comparison of the levels thus estimated and known levels at the four standard ports in the study region (Immingham, Cromer, Lowestoft, and Felixstowe) suggests that this approximation tends to slightly overestimate (order of 10 cm) the level of HAT. The dependency of intertidal wetland development on tidal range is well documented, with equivalent landforms occurring at higher levels relative to the mean sea level in macrotidal settings than in microtidal ones [27,36]. The elevation range available for the development of intertidal landforms can be described as an 'accommodation range'. The DEM elevations above 0 m ODN were therefore normalised to reflect the spatial variability in accommodation range, which was itself approximated as 0.5× MSTR. This produced a normalised elevation (NE) raster that did not exhibit substantial dependency on tidal range and was used to scale the resulting cost functions for the intertidal zone.
A 'supratidal' zone was defined as areas above HAT but seaward of the SMP2 vector marking the terrestrial limit. Such areas may still be inundated and therefore permit water flow during exceptional events, so this zone was given a uniform but very high cell cost relative to the other zones. A 'terrestrial' zone landward of engineered defences was modelled as impassable for tidal waters.
The zones were combined to produce an overall cost surface raster which was converted to a cost-distance raster where each cell value represents the cumulative cost to reach it by the least-cost path from the 10 km offshore limit of the domain. The weightings for each zone were manually adjusted to produce the best achievable visual replication of expected flow routes for tidal waters to reach a given location through complex channel networks. Figure 3 shows a schematic cross-shore transect from the engineered defence to the 10 km offshore limit of the domain (top panel). Per-cell and cumulative costs are plotted in red, with the levels of 0 m ODN and HAT also indicated. The four zones upon which the overall cost surface was based are indicated along with their final weighting functions (where D is the horizontal distance in metres and NE is the normalised elevation). The lower panel of Figure 3 shows examples of the least-cost paths calculated for random pixels superimposed on both the local per-pixel cost surface and aerial photography to demonstrate the ability of the method to reproduce expected flow routes.

Distance from Creek/Edge of Marsh Parcel
This measure represents the well-established tendency for sediments to be deposited rapidly once creek banks become overtopped during inundation, leading to the development of creek levees and limited deposition rates in areas away from cliffed seaward margins [37,38]. To represent this dynamic, a measure of distance from a creek or the edge of a marsh parcel was used. To facilitate analysis at scales commensurate with the size of Landsat pixels used for other variables (30 m by 30 m), a reduced resolution product representing offshore or large-creek margins of marsh parcels was derived from the Environment Agency's salt marsh extent layer [17]. Large creeks were defined as those wider than the diagonal of a Landsat pixel (ca. 42.5 m), since these are the units for which internal morphological changes were measured. Creeks and pools narrower than 42.5 m were removed using a 21.25 m outer buffer on the polygons of the salt marsh extents layer, followed by the dissolving of any overlaps thus created and 21.25 m inner buffer operations. This spatial resolution limit, imposed by the buffering process, effectively disregards the smaller creeks, which nevertheless perform important system functions. This fact must be taken into account when interpreting the findings, and we comment further on this aspect in the discussion. The other metrics used, including δP veg , are not sensitive to spatial resolution limits in the same manner.
Landward limits of the marsh polygons were removed by extracting only those portions of perimeter lines that were greater than 30 m (1 landsat pixel) from the SMP2 vector, ensuring that distances do not represent those from the landward limit. For each parcel, the euclidean distance from the edge was then calculated on a 30 m grid aligned to that used for δP veg estimation. The metric is henceforth referred to as Euclidean Distance (ED).

Integrity of Marsh Platform
The first morphological metric describes the integrity of the marsh platform by the proportion of each 30 m satellite pixel's area that contains undissected vegetated marsh areas. It is referred to as pCore.
The Morphological Spatial Pattern Analysis (MSPA- [39]) tool, supplied within the GUIDOS toolbox [40], is a morphological segmentation algorithm designed to assess habitat connectivity and fragmentation-for example, in forest ecosystems [41,42] or for planning of green infrastructure [43]. MSPA operates on a binary raster of 'foreground' and 'background' pixels, where foreground represents the habitat or land cover type of interest. It classifies foreground pixels according to how many other foreground pixels they are adjacent to and a distance parameter controlling the width of areas considered to be 'egde' because they are close to background pixels. Pixels are allocated to one of seven foreground classes (core, edge, perforation, bridge, loop, branch, and islet) and a background class. Each of these classes can also have an attribute describing whether it is entirely surrounded by other foreground classes (internal) or with connectivity through adjacent background classes to the edge of the raster (external). Loops and bridges can additionally be defined as appearing independently or within edges or perforations. The result is 22 possible types of feature that are extracted from the binary raster. In this context, foreground pixels represent vegetated marsh platform areas while background pixels correspond to bare sediments, channels, and pools. Although these classes were originally designed to describe habitat connectivity features, analogues can be readily defined for most classes in the context of fine-scale saltmarsh morphology (Table 1). MSPA was applied to a 25 cm-resolution binary raster derived from the Environment Agency's salt marsh extents layer [17], which describes recent (2008-2010) marsh morphology. Core areas were defined as parts of vegetated marsh greater than the MSPA edge width parameter (2.5 m) from a background area. The pCore metric was the aggregation of both internal and external MSPA 'core' classes for a given satellite pixel area.

. Limitations to Tidal Connectivity
The second morphological metric describes the proportion of pixel area comprised of unvegetated areas and their periphery that are unconnected to the drainage network; it is denoted pUncon. It is also derived from the MSPA segmentation and reflects pools or pans and their periphery that do not have connectivity via channels to the offshore domain. The pUncon metric was computed as the aggregation of the 'internal' MSPA classes, including unvegetated areas.
An example area of the MSPA segmentation for six adjacent satellite pixels is shown in Figure 4, with the respective values for the pCore and pUncon metrics. These two morphological metrics are interpreted as describing the marsh morphology within a dataspace that separates, at its extremes, between contiguous, undissected marsh platform (high pCore), marsh areas that are heavily fragmented by creek networks (low pCore, low pUncon), and areas that have a high density of salt pans and pools (low pCore, high pUncon). The inset graph on Figure 4 illustrates the position of the depicted pixels within this data-space. The processes leading to these morphologies are likely to be different, as are their likely future evolutionary trajectories. Creek networks, particularly low-order ones, are typically stable in the platform once established [44], and maintain a dynamic equlibrium with the tidal prism [45]. Interior areas of intact marsh platforms may remain so or become punctuated by pools or pans [46][47][48], with the density of pans being related to the tidal range and sediment type at a national scale in the UK [49] and the elevation and distance from marsh edge at a site scale [47].

Probability of Observing Marsh Degradation
Marsh degradation was considered to have occurred for pixels where the δP veg was negative. For each of the topological and morphological metrics, the probability of a location exhibiting degradation was computed from the proportion of pixels where δP veg < 0 in each of 100 equal intervals across the observed ranges of the four topological or morphological metrics described above (CD, ED, pCore, pUncon). Least-squares regression analyses were used to establish relationships for each of the metrics between the observed probabilities of degradation and interval centre points.

Morphodynamic Change
δP veg was estimated for approximately 180,000 30 m pixels, amounting to approximately 16,200 ha of marsh. An example of the resulting metric for about 600 ha of salt marsh in Hamford Water, Essex, is presented in Figure 5 to illustrate the degree of spatial variability observed in changes to marsh integrity over relatively small spatial scales. The inset histogram shows the distribution of estimates across the entire study domain from the Humber to the Thames, which can be seen to be normal but with a mean somewhat above zero, suggesting a dominance of vegetation establishment within the domain, even though this is not necessarily evident within the Hamford Water subset depicted. The overall probability of degradation for the entire dataset was 0.144 for the approximate time period 1985-2016.

Topological and Morphological Relationships
No significant regression relationship was found between Cost Distance (CD) and the probability of degradation when using the entire dataset. The CD parameter, however, has very long tails to its distribution, particularly at the upper end of the scale. At the lower end, these pixels represent the seaward limit of open-coast marsh areas, while at the upper end they are mainly high, supratidal regions towards the heads of estuaries or around islands of high ground. These locations may not be expected to experience the same controls on marsh morphodynamics as the majority of interior and estuarine marsh areas, so might confound the regression analysis for such areas. Areas with exceptionally low CD values coincide with the seaward margins of open coast marshes, which have the potential to retreat [14] and advance rapidly [24]. Values of δP veg in such areas therefore tend towards extreme values and reflect changes in marginal position more than changes in marsh platform integrity. Where CD is exceptionally high, δP veg values are associated with the upland boundary and may again not reflect changes in marsh platform integrity but rather shifts between terrestrial and halophytic vegetation communities, which may have completely different reflectances. The extremities of the CD distribution were discarded where CD ≤ 500,000 and CD ≥ 1,400,000. 91% of the population was retained, and all the locations excluded within Hamford Water are marked, for illustration purposes, as black pixels within the red circles on Figure 5. With the extremities of the distribution excluded, a significant positive linear relationship between CD and the probability of degradation was found (Adjusted R 2 = 0.37, p ≤ 0.001-see Figure 6, Panel A). The low positive coefficient suggests that, as CD increases-implying a reduction in tidal connectivity to offshore waters-the probability of marsh degradation increases slightly. CD values do not represent real-world physical quantities. For instance, a doubling of CD could occur because of multiple different combinations of increased horizontal distance and/or elevation. The precise coefficients of this relationship are therefore not readily interpretable beyond their sign and approximate magnitude.
A second-order polynomial relationship was found between Euclidean Distance (ED) and the probability of degradation, whereby areas with a low ED exhibit moderate probabilities of degradation in the order of 0.1, with this probability declining to a minimum of less than 0.05 at around 350 m before increasing again with ED to values in excess of 0.3 at distances of 1200 m. The relationship is significant at p ≤ 0.001, and the high adjusted R 2 of 0.76 suggests that ED is an important determinant of changes in condition ( Figure 6, Panel B). The ED metric is directly comparable to that used by Kearney and Rogers [20] to predict changes in marsh condition, and the regression fit from their study is superimposed on the figure for comparison.
Significant relationships were also found between both pCore and pUncon and the probability of degradation ( Figure 6, panels C and D, respectively). The former relationship implies a decrease in the probability of degradation with increasing values of the pCore, although very high values of pCore become associated with increased probabilities of degradation again. With high adjusted R 2 values of 0.63 for pCore, the proportion of a pixel that is core vegetated marsh area is strongly related to changes in marsh condition. Pixels containing salt pan features were relatively rare within the dataset, with n = 13157, representing 9.2% of all pixels. The regression only included those pixels where pUncon > 0. The relationship is significant with an R 2 value of 0.63, suggesting a strong positive relationship between the proportion of an area that comprises salt pans and their peripheral marsh features and the history of morphological evolution.
Summaries of all four regression relationships are provided in Table 2.

Morphodynamic Change
The regional data for the Humber to Thames study area show an overall trend towards increasing marsh platform integrity over time, with rapid increases in integrity being much more common than rapid losses (inset to Figure 5). This pattern is, in part, a result of the rapid infilling and marsh advance that occurred in The Wash embayment during the observation period , where rates of margin advance up to 75 m per year were observed between 1992 and 2014 [24]. This established a large area of marsh that is included in the analysis and is associated with very high values of δP veg . The pattern of increase and decrease in internal platform integrity is much patchier elsewhere, as seen for Hamford Water ( Figure 5).
A degree of bias towards positive values of δP veg was introduced by the choice of a marsh presence/absence mask from relatively late in the study period. Phelan et al. [17] base their mapping on data from 2008 to 2010. This implies that areas that may have contained marsh early in the study period, but which became completely unvegetated prior to 2008 and would produce negative δP veg , will have been excluded from the analysis. A visual assessment of the aerial photography from 2013/14 and 1992 suggests that this phenomenon was rare within the region. The scale of negative changes in the overall horizontal extent of the marshes over the study period, relative to the Landsat pixel size, is typically small. Only isolated pixels at the seaward limits of some open-coast marshes are likely to have been affected by this bias [24]. Our findings support the idea that vegetated areas of marsh platform have the potential to increase in platform integrity much more rapidly than they deteriorate. This is likely the result of the potential for the rapid colonisation of large horizontal extents of unvegetated surfaces that accrete to a critical elevation to support seedling survival [29] near-simultaneously because of their very low gradients. Conversely, during erosive phases sediment loss is likely to be more gradual due to biostabilisation [50,51], and is also likely to be more localised in areas of high hydrodynamic stress. Assuming that morphological change can be inferred from vegetation extent change, the data suggest that, in marsh interiors at least, accretionary processes tend to outpace erosional ones.
We address change over a single time period only. Our methodology is not, therefore, able to detect regime shifts or abrupt changes as distinct from more gradual trends that result in the same δP veg value. A natural development of our work, as the data archives become longer and denser, would be to incorporate temporal segmentation approaches similar to those used by the LandTrendr algorithm [52,53]. This would allow us to investigate whether the topological and morphological controls on abrupt regime shifts differ from those on gradual changes. Such work is beyond the scope of the current study.

Geomorphic Setting (CD)
No significant relationship was found when considering probability of degradation across the entire range of CD. This is likely because CD has a large range between 3.9 × 10 5 and 3.1 × 10 6 , with over 90% of values falling in the interval 5.0 × 10 5 and 1.4 × 10 6 . The upper end of the distribution in particular is therefore very sparse. This leads to very small numbers of pixels being represented in each of the 100 intervals across the CD distribution. With the relatively low prevalence of degradation in the order of 10% within the entire dataset, these small samples at high CD values often fail to represent any degraded pixels, resulting in probabilities of degradation of zero. Even if some degraded pixels are sampled, the resulting probability estimate is unlikely to be usefully precise while the sample size remains below about 50 pixels. Furthermore, the areas where CD> 1.4 × 10 6 are typically high-elevation areas which are likely to be supratidal and may therefore be expected to be responding to different processes controlling their dynamics when compared to the intertidal marsh areas that are of principal interest here. Similarly,the data density is very low where CD < 5.0 × 10 5 , as these areas represent the seaward limits of open-coast marshes. The observations in these areas are likely to be dominated by the wave-driven erosion of the marsh margin rather than the internal marsh changes that δP veg represents. The probabilities of degradation in these areas are relatively high, aligning well with the observations of the widespread retreat of open coast marsh margins in the region [24]. Excluding such areas from the analysis is therefore justifiable, as it enables insight into dynamics in the intertidal zone that are otherwise masked by variance in other areas that are either artefactual or likely to arise from a different suite of processes to those types of marsh interior processes that this study seeks to assess.
When the range of CD was limited to the interval 5.0 × 10 5 and 1.4 × 10 6 , the regression analysis showed a significant positive linear relationship between CD and the probability of a pixel showing degradation. Low CD values are associated with minimum probabilities of degradation (Panel A, Figure 6). This implies that the interiors of open-coast marshes experience little degradation, while those towards the head of estuaries or tidal inlets are more susceptible. This aligns well with the conceptualisation of sediment sources for this study region, which is dominated by minerogenic marshes [26] and offshore tidal waters are the primary source of allochthonous inputto support marsh elevations [54]. Those areas of marsh that are least tidally connected (high CD) are shown to be the most vulnerable, since the delivery of sediments is most restricted in these areas. This finding contrasts with those of Kearney and Rogers [20], who observe an increasing probability of degradation at shorter distances up-estuary. This is because the domain that they studied on the east coast of the US has much more significant riverine inputs and the locations of sedimentation are therefore fluvially, rather than tidally, dominated. Additionally, the Chesapeake Bay (US) marshes are highly organogenic, with the upper marsh sediment total organic content being reported to be around 80% at Rhode River [55], which implies that autochthonous accumulation dominates. By comparison, the UK east coast marshes rely more heavily on allochthonous inputs with much lower total organic contents than those of the US, at around 15% [56]. The relationship we find between CD and the probability of degradation is relatively weak (R 2 = 0.37), reflecting the fact that other factors, likely occurring at smaller scales, may be more significant. These may include the controls represented by the other ED, pCore, and pUncon metrics in this study, but also factors such as proximity and the magnitude of fluvial sediment loads that we have been unable to address here. The strength of the relationship may also arise, in part, from the failures of the simplified CD metric to perfectly represent tidal connectivity within the domain. Nevertheless, given the size of the region considered, the heterogeneity of marsh settings represented and the methodological limitations, the fact that any relationship emerges implies that the position of a marsh within the wider coastal setting is an important control on its vulnerability to sediment starvation and subsequent degradation. The contrast between our findings and those of Kearney and Rogers [20] suggests that the marsh position is not, however, a diagnostic parameter, since its effect varies between regions. Rather, it modulates the impacts of larger-scale geomorphological contexts (that may be considered as boundary conditions at the scale of the analysis presented here), such as sediment sources and delivery pathways. This finding implies that empirical attempts to predict morphological change in coastal wetlands must be nested within a hierarchy of process understanding in order to ensure that outcomes are appropriate for a particular location. Not all regions behave equally.

Distance from Creek/Edge of Marsh Parcel (ED)
The relationship established between Euclidean Distance (ED) and probability of degradation is significant and stronger (R 2 = 0.76) than that found for CD. It is also nonlinear with initially high probabilities at short distances that decline to a minimum around 350 m before increasing dramatically. There is a notable similarity with the relationship observed in Chesapeake and Delaware Bays (black dashed line on Figure 6). The initial decline observed here is somewhat less rapid than that found by Kearney and Rogers [20], and beyond 400 m there is a strong increase in the probability of degradation with increasing ED at a rate that exceeds that found by Kearney and Rogers [20].
The inter-regional similarity observable regarding ED suggests that topological factors controlling marsh stability at scales below those represented by CD (i.e., within-marsh scales rather than within-estuary scales) may be operating in very similar fashions between regions, although the differences in the rates of decline and increases in degradation probability around the 350 m minima suggest a possible dependence on the scale/size of the marsh systems being observed. Marsh parcels are substantially smaller in the UK than in the USA example, where ED values of up to about 100 km were observed, although the goodness of fit above about 1 km declines. Differences in drainage network morphologies may also contribute. For the context of the east coast of the UK, the initial decline in probabilities can be interpreted, in part, as a function of the scale of channel identified in the basemap used for the distance calculations. Only major creeks exceeding 42.5 m width were included. Where the ED values are very low, the degradation measured by δP veg could have a component that reflects the marginal erosion of the banks of large channels, where fetch is sufficient for wave-driven erosion to become significant in determining the spectral change within the pixel over the study period. Thereafter, however, it still appears that marsh vulnerability decreases with increasing distance from these major channels. This could be associated with the increasing drainage density of smaller creeks (not resolved by the ED metric used here), leading to shorter unchannelised lengths, a greater preponderance of levee-effects [19], and ultimately a greater potential for the drainage of and allochthonous sediment import to marsh areas within 350 m of a major channel. Beyond this distance, the probability of degradation increases again, which we attribute to declining drainage density and efficiency, leading to decreasing potential for tidal flushing, drainage, and sediment import as distances from the primary tidal channels continue to increase. Internal marsh areas without efficient drainage and a long way from major channels tend to experience relatively little elevation gain from external sources [57] and, particularly in the context of relative sea level rise, vegetation experiences increasing water and salt stress, potentially leading to die-back and substrate collapse [58]. This phenomenon is observable within the region and is illustrated by the hotspots of negative δP veg in marsh interior areas shown in Figure A1 within the Supplementary Material.

Percentage Core Areas (pCore)
pCore is based on a 25 cm spatial resolution and therefore resolves differences in marsh morphology at a sub-pixel scale in terms of δP veg . The same is true for pUncon. The relationship observed whereby the probability of degradation decreases to a minimum at pCore values of around 70% before levelling off or increasing slightly implies that fragmented areas of marsh (low pCore) reflect a history of greater marsh degradation. This implies that more degraded areas of marsh may continue to degrade at a faster rate than areas of platform with a high integrity. As a marsh platform becomes fragmented, pCore declines and the perimeter-area ratio of the remaining marsh parcels increases, providing a greater effective marginal length that may be vulnerable to erosion by hydrodynamic forces. The stability of the morphological state therefore decreases with increasing fragmentation, representing a positive morphodynamic feedback controlling marsh degradation. The slight increase in probability at very high pCore values may reflect the tendency for marsh areas that are separated from minor creeks by only a few metres to begin to degrade as a result of water stress [58], which Ursino et al. [59] showed can increase rapidly away from channels. Additional contributors to marsh degradation occurring at high pCore values may be reduced sediment inputs [19,57] or increased salt stress [60] in the interiors of marsh parcels.

Percentage of Unconnected Areas (pUncon)
The relationship between pUncon and the probability of degradation is significant, with a R 2 value of 0.63. The form of the relationship takes an almost opposite form to that for pCore, with higher values of pUncon being associated with higher probabilities of degradation. pCore and pUncon are somewhat negatively correlated to each other (R = −0.11, p ≤ 0.001), so this inverse relationship may in part reflect the same processes as are discussed with reference to pCore. The analysis for pUncon, however, incorporates only those marsh pixels where pans or pools are present, which represents a small subset of the dataset (9.2%). As such, this analysis allows for the isolation of behaviours exclusively in areas where pans are observed. Areas with a higher pUncon and therefore a greater proportion of pan-related features exhibit higher probabilities of degradation. Notably, for fragmentation identified by the pUncon metric, the maximum probabilities of degradation observed are much higher (≈0.5) than those calculated on the larger dataset of pCore (≈0.2), where low metric values also indicate fragmentation. Low pCore values, by contrast, largely reflect fragmentation by creek networks and bare surfaces connected to them. The difference in magnitudes of probability identified by the two regression relationships implies that the presence of pans indicates a much more vulnerable landscape configuration than one fragmented by creeks. Neumeier et al. [61] identified the life cycle of pans in Canadian marshes whereby they undergo a phase of active expansion before achieving maturity. The data presented here also support the conclusion that pans have a tendency to expand over decadal timescales. The analysis cannot identify causality, however. It is possible that the presence of pans establishes a feedback encouraging their subsequent expansion as the marsh develops [48]. Questions also remain over the causes of initial pan formation, with some theories suggesting that event-based phenomena, such as smothering of vegetation by rafted debris, may cause pans to develop. French et al. [62] noted that pan densities are lower on the more enclosed of the backbarrier marshes in Norfolk, where less debris may be rafted into the marshes over shorter seaward margins. Their finding is correlative but Pethick [47] argues that a mechanism for pan initiation on mature marshes must exist. Whatever the cause of pan initiation, our findings suggest that, once established, marsh areas containing pans are highly vulnerable compared to other areas. It is possible, therefore, that stochastic events, such as debris rafting, could initiate a long-term landscape vulnerability. In contrast, pans may simply be diagnostic of other conditions such as relative sea level changes, tidal range, or coastal configuration, rather than morphodynamic drivers in their own right [49]. If pans are not, in themselves, morphodynamic drivers, then our findings would suggest that the external controls causing pan formation continue to cause subsequent expansion and the conclusion that panned landscapes are vulnerable stands. Further work is needed to address the question of controls on salt pan development to better understand the fundamental causes of this vulnerability.

Conclusions
Our findings show that there exist certain overarching controls on the vulnerability of salt marshes to degradation at a range of scales. We demonstrate this observation through the application of a systematic spatial analysis of a simple, satellite-image-derived measure of change. By applying this methodology at the regional scale and across a range of estuarine, backbarrier, embayment, and open coast settings, we illustrate the nature of these controls more clearly than has hitherto been possible.
This work has presented an inherently scalable method for monitoring the platform integrity of salt marsh surfaces at a sub-pixel scale that is easily adapted to incorporate current and future data sources (such as imagery from the more recently launched Sentinel 2 satellites). The δP veg metric presented here is a continuous variable representing change in marsh surface integrity over a given time period. It therefore provides substantially more statistical possibilities than the MSCI [20] or than we have explored in this study.
We express topological and morphological factors at a range of scales as metrics describing the connectivity of tidal waters, position within marsh, marsh platform integrity, and the prevalence of salt pans or pools. We demonstrate that all of these are significantly related to the morphodynamic evolution of marshes. This spatial hierarchy of morphological controls may continue to be applicable at both larger and smaller scales than are addressed here. The presence of controls at larger scales is shown by the finding of contrasting relationships between geomorphic setting and marsh degradation in this study and that of [20]. We hypothesise that this difference arises from regional-scale contrasts in climate, biota, sediment supply, and sea-level history that cause fundamental differences in the processes by which the marsh systems function and evolve. At a whole-marsh scale, however, responses to the distance from the marsh margin seem to be fairly consistent between regions, suggesting commonality in some of the processes governing marsh dynamics. At a smaller scale than is considered here, a variety of factors may become important, such as the spatial distribution of vegetation types [46,63] or phenotypes [64], which may determine sediment erodibility.
We have presented aggregated relationships that emerge at a regional scale. Overall, these suggest that marsh areas that are already exhibiting some form of fragmentation that lie far from the nearest creek and towards the heads of estuaries and inlets are the most likely to exhibit degradation. The prediction of changes in a particular location remains challenging. No relationships were found between the continuous values of δP veg itself and any of the predictors presented here using basic statistical techniques. Only when the dimensionality of the dataset is reduced to probabilities do aggregated relationships emerge. The large number of parameters influencing marsh dynamics, the multiple scales of processes involved, and the complex interactions between them means that statistical prediction of internal marsh changes will require methods capable of capturing the interactions and high dimensionality that are inherent to such natural systems. Further work is needed to explore the potential of more sophisticated statistical techniques drawn from machine learning domains to synthesise datasets such as the one presented here into location-specific values of parameters describing the sign and magnitude of predicted marsh evolution.

Abbreviations
The following abbreviations are used in this manuscript: The Normalised Difference Vegetation Index (NDVI), a dimensionless index taking values between −1 and 1, is a proxy for the amount of chlorophyll (and therefore vegetation) present, with higher values indicating more chlorophyll. The NDVI is the normalised ratio between reflectance in the red and near-infrared bands. Healthy vegetation produces positive index values (approaching one) [65]. Conceptually, the NDVI for a pixel in a satellite image can be considered to be a function of the percentage of that pixel covered by vegetation and the nature of the vegetation within that pixel. Applying an assumption of approximate stationarity in the reflectance of the bare sediments, changes in the percentage vegetation cover within the pixel are reflected in a change in NDVI, with the signal potentially being modulated by any attendant changes in the community composition (and therefore spectral signature) and vigour of the vegetation present.
Summer NDVI is relatively insensitive to the successional stage of the vegetation within this study region, while being lower for pioneer stages than mature stages during winter. Field spectroscopy was conducted at Tillingham, Essex, throughout the period 2015-2016. A total of 58 pioneer spectra and 103 mature marsh spectra were collected using an SVC HR1024i spectroradiometer and calibrated reference panel. The resulting hemispherical-conical reflectance factors were convolved to Landsat-8 band responses, from which NDVI was calculated. Pioneer vegetation (dominated by Salicornia and Spartina species) has similar Landsat-8 NDVI in summer (Jun-Oct) to the perennial marsh canopy (t-test, p ≥ 0.05). However, the near-complete dieback of pioneer vegetation leads to the exposure of bare sediment and therefore lower winter NDVIs.
In Google Earth Engine [32], all scenes intersecting the study area were extracted from the Landsat archives for Landsat-5, -7, and -8 at a spatial resolution of 30 m, providing coverage between 1985 and 2018 (excluding Landsat-7 images subsequent to the failure of the scan line corrector on 31 May 2003). Level 1 Top-of-Atmosphere images were imported and filtered for cloud cover based on their metadata cloud score and subsequently on a per-pixel basis using the Simple Cloud Score algorithm. Thresholds for excluding data were set at 20% for both stages. NDVI was computed for all the remaining pixels of all images.
Differences in sensor response, and the NDVI values derived from them, are widely recognised [66] and have been shown to produce bias in time series analysis if combined without corrections being applied [67]. NDVI images were corrected for different sensor responses to produce a consistent time series of Landsat-7-equivalent NDVI values. Landsat-8 scenes were corrected using the cross-calibration derived by Roy et al. [68]. No existing cross-calibration was available for Landsat 5. To calculate a calibration function pixels were identified that fell within large areas of established marsh (those exceeding 1500 m 2 in the UK salt marsh extent map [17]) and more than 42.4 m (the diagonal of a 30 m Landsat pixel) from the edge of the marsh area. These were considered likely to be 'pure', fully vegetated locations. All the Landsat 5 and Landsat-7 scenes were selected for the period between the launch of Landsat-7 on 15 April 1999 and 31 May 2003. This process resulted in a total of 3747 pixels being selected for the cross-calibration using 5 Landsat-5 scenes and 40 Landsat-7 scenes. The per-pixel mean NDVI values were compared using a second-order polynomial, since a degree of saturation occurred in the Landsat-5 sensor when compared to Landsat-7 (p ≤ 0.0001, R 2 = 0.756). The cross-calibration functions applied are detailed in Table A1. Table A1. Cross-calibration functions to correct NDVI from Landsat-5 (L 5 ) and Landsat-8 (L 8 ) [68] to Landsat-7 (L 7 ) equivalence prior to time series analysis.

Transformation
Function R 2 L 5 to L 7 L 7 = 0.085 + (0.460L 5 ) + (0.247L 2 5 ) 0.756 L 8 to L 7 L 7 = −0.011 + 0.969L 8 0.906 The corrected Landsat-7-equivalent time series of NDVI images was processed in the Google Earth Engine to produce per-pixel linear trend coefficients. These trend coefficients were weakly related to the mean NDVI for the pixel over the entire time series (R 2 = 0.22, p ≤ 0.001). This was interpreted as arising either from a CO 2 fertilisation signal causing vegetation in general to become 'greener' throughout the time series. Alternatively, it may reflect a pixel-scale morphological control on change in vegetation distribution, whereby pixels with higher platform integrity tend to increase their integrity more than those with lower integrity. Either of these effects would confound the investigation of topological controls. The relationship was therefore removed by using the residuals of the linear fit in Equation (A1) where δN/δt denotes the NDVI trend andN denotes the mean NDVI. The residuals were subsequently standardised.
δN/δt = 0.011254N − 0.001430 (A1) Standardised residuals were converted into sub-pixel estimates of percentage change in vegetation cover (henceforth δP veg ) by calibration against geo-referenced vertical aerial photography supplied by the UK Environment Agency. The observed range of the standardised residuals was divided into 100 equal intervals and one pixel was selected at random within each interval. Some intervals towards the tails of the distribution contained no pixels, precluding their inclusion. In total, 83 pixels were selected. Pixels were clipped from the earliest available aerial photography (1992, panchromatic, 25 cm resolution) and from colour photography from 2013/2014 (20 cm resolution).
In Matlab, 50 random points were generated per Landsat pixel image pair (1992 and 2013/14), and sequentially superimposed on each image. A single operator visually assessed the surface cover and allocated each point as either vegetated or not vegetated. Manual assessment and attribution is assumed to provide the highest achievable accuracy in this context. A total of 4150 points were manually attributed for each year. From the proportions of vegetated and unvegetated points in an individual Landsat pixel area, the percentage vegetation cover was estimated. The difference between the percentage of the pixel area that was vegetated in 2013/14 compared to 1992 was calculated. A linear relationship was found between the standardised residuals and the change in vegetation cover. The linear relationship between the change in percentage vegetation cover between 1992 and 2013/14 estimated from aerial photography and the standardised residuals of the NDVI trend analysis (R trend ) is given in Equation (A2). The R-squared value of 0.87 implies that the change in vegetated area is the major contributor to changes in the residual trend. The unexplained variance is likely the result of other factors not accounted for by vegetation extent alone, such as species compositional change.
The method outlined above produces estimates of change in marsh platform integrity as continuous values between ±100%. This offers many statistical possibilities that will be explored in future work. For the purposes of the current study, estimates were thresholded at zero to produce a binary gain/loss indicator allowing for the alignment of statistical methods with the work of Kearney and Rogers [20], who discerned between pixels that were either degraded or not degraded. The performance of the calibrated δP veg method was evaluated at a more local scale than the entire study domain in case large scale spatial dependency associated with, for example, latitude was a significant factor affecting calibration. The area selected for validation was Hamford Water, Essex, a tidal inlet containing approximately 600 ha of marsh that has been shown to exhibit a wide variety of morphodynamic behaviours in close spatial proximity to each other [69]. Analysis up to the point of manual calibration was repeated using only Hamford Water as the domain to establish locally derived residual trends which were classified into ten classes using a Jenks Natural Breaks algorithm. The largest contiguous area for each class was identified and from that area ten pixels were selected at random for validation, providing 100 pixels representing a range of morphological behaviours. These were analysed by manual point attribution from aerial photography following the same method as previously outlined. The change in percentage