Identification of typical ecohydrological behaviours using InSAR allows landscape-scale mapping of peatland condition

Better tools for rapid and reliable assessment of global peatland extent and condition are urgently needed to support action to prevent further decline of peatlands. Peatland surface motion is a response to changes in the water and gas content of a peat body regulated by the ecology and hydrology of a peatland system. Surface motion is therefore a sensitive measure of ecohydrological condition but has traditionally been impossible to measure at the landscape scale. Here we examine the potential of surface motion metrics derived from satellite interferometric synthetic aperture radar (InSAR) to map peatland condition in a blanket bog landscape. We show that the timing of maximum seasonal swelling of the peat is characterised by a bimodal distribution. The first maximum, usually in autumn, is typical of “stiffer” peat associated with steeper topographic gradients, peatland margins, and degraded peatland and more often associated with “shrub”-dominated vegetation communities. The second maximum, usually in winter, is typically associated with “softer” peat typically found in low topographic gradients often featuring pool systems and Sphagnum-dominated vegetation communities. Specific conditions of “soft” and “stiff” peats are also determined by the amplitude of swelling and multi-annual average motion. Peatland restoration currently follows a re-wetting strategy; however, our approach highlights that landscape setting appears to determine the optimal endpoint for restoration. Aligning the expectation for restoration outcomes with landscape setting might optimise peatland stability and carbon storage. Importantly, deployment of this approach, based on surface motion dynamics, could support peatland mapping and management on a global scale.

1 Introduction

The conservation of well-functioning peatlands and restoration of degraded peatlands to reduce and ultimately mitigate land-use-related emissions of atmospheric carbon dioxide is now a global priority (Leifeld and Menichetti, 2018; Amelung et al., 2020; Günther et al., 2020). Furthermore, to support the implementation of national peatland management plans and restoration initiatives, cost-effective measures to record current peatland condition and restoration progress are urgently required (Crump, 2017). Mapping peatland extent and condition has long been recognised as a huge challenge over large, remote, wet, and often discontinuous peat-forming regions where field-based surveys are impractical and expensive (Lees et al., 2018). Alternatives such as thematic mapping based on optical remote sensing (visible and near-infrared) are increasingly used (Artz et al., 2019; Minasny et al., 2019; Lees et al., 2020), but the number of observations in regions with frequent cloud cover such as peatland in the northern latitudes and the tropics reduces the number of possible surface observations. Radio detection and ranging (radar), which is sensitive to physical properties of the surface, provide an effective, more frequent option, given that microwave frequencies can penetrate cloud cover and return a measured signal from the ground (Minasny et al., 2019; Poggio and Gimona, 2014). For example, using the ESA Sentinel-1 synthetic aperture radar (SAR) satellites it is now possible to observe a peatland surface anywhere at high frequency (6 to 12 d) with continuous spatial coverage. When this is combined with the technique of SAR Interferometry (InSAR), it allows detection of surface displacement, an indication of peatland condition, as a time series of observations (Sowter et al., 2013).

In peatland, the rise and fall of the surface, sometimes described as “bog breathing” (Kulczynski, 1949; Baden and Eggelsmann, 1964; Mustonen and Suena, 1971; Hutchinson, 1980; Kurimo, 1983; Almendinger et al., 1986; Price, 2003; Price and Schlotzhauer, 1999), is one of the key self-regulating feedback mechanisms in peatland, providing resilience and maintaining function during periods of hydrological stress (Money and Wheeler, 1999; Waddington et al., 2015; Mahdiyasa et al., 2022). This “surface motion”, which is a poro-elastic mechanical response to ecohydrological processes, results from the collapse and expansion of large pores in response to changes in the mass of water stored and associated stresses within the peat (Price, 2003; Mahdiyasa et al., 2022). Mechanical deformation of the peat body and consequent surface motion can also modify the ecohydrology of a peatland via compaction, slope failure, and pipe formation (Waddington et al., 2010, 2015). Small-scale field observations indicate that peat surface motion is influenced by changes in water level (Roulet, 1991; Price, 2003; Kennedy and Price, 2005; Fritz at al., 2008; Alshammari et al., 2020), vegetation composition (Howie and Hebda, 2018; Alshammari et al., 2020), micro-topography (Waddington et al., 2010), accumulation and upward migration of methane bubbles (Glaser et al., 2004; Reeve et al., 2013), and land management (Kennedy and Price, 2005).

Collectively these results suggest that peatland surface motion could be a sensitive indicator of peatland function on a landscape scale. So far, InSAR investigations have focused on discrete, small-scale (<1 km²) peatlands (Fiaschi et al., 2019; Tampuu et al., 2020), identifying the potential range in timing and amplitude of seasonal peatland surface motion (Alshammari et al., 2018, 2020) and its relationship to precipitation (Fiaschi et al., 2019), water level (Alshammari et al., 2020; Tampuu et al., 2020), and vegetation composition (Alshammari et al., 2020). However, peatland landscapes contain a continuum of topographic, ecological, hydrological, and management regimes, and these small-scale studies have not captured the full spectrum of peatland conditions between degraded and near natural.

In this paper, we determine whether surface motion measured by InSAR can be used to quantify peatland condition continuously over a complex peatland landscape. Using the APSIS (Advanced Pixel System using Intermittent Small Baseline Subset, formerly known as the Intermittent Small Baseline Subset (ISBAS)) method, which is capable of generating spatially continuous measures of vertical surface motion over peatland (Sowter et al., 2013), we measure time series of surface motion over our study site at a high spatial and temporal resolution. Specific time series metrics are then compared to independent measures of peatland condition to determine their relationship. By doing this we relate surface motion metrics to the continuum of ecohydrological conditions in this peatland landscape. Finally, we demonstrate how surface motion metrics can be used to map the ecohydrology of a peatland system. By doing so we illustrate how our new approach could be applied to monitoring the response of global peatland to restoration, management, and climate change.

2 Study location

The Flow Country peatlands, northern Scotland (Andersen et al., 2018), exist in a range of topographic, hydrological, and management settings, leading to a range of different conditions, e.g. highly eroded uplands to relatively intact low-lying peatland with pool systems, with activities such as forestry, drainage, and grazing superimposed. Our chosen study site is 930 km² of blanket bog, ranging from 50 to 600 m a.s.l. (Fig. 1). From the 19th century onwards, management has involved burning to support grouse shooting, artificial drainage of the driest areas of peatland, subsidised agricultural improvement, and later afforestation programmes (Sloan et al., 2018). More recently, near-natural areas have been designated for conservation (Lindsay et al., 1988), and previously afforested and drained areas started undergoing restoration. Some areas are also actively eroding, particularly at the highest altitudes (Hancock et al., 2018). This complex mosaic of near-natural and modified peatlands makes the study site particularly suited to an investigation on the use of InSAR for mapping the peatland condition. To understand the variation and distribution in the characteristics of the APSIS-derived time series, we analysed five well-documented 10 to 15 km² sub-sites within the study area (Table 1; Fig. 1).

Figure 1 The study location (inset) and true colour satellite image composite covering the study area, outlined, in the Flow Country, northern Scotland. Forested areas are dark green; more intensive agriculture appears as lighter greens to the north-east, with the remainder of the image mostly consisting of blanket bog. Superimposed on the image are the main river network and the location of pool features shown with sub-sites for detailed analysis marked as SS1 to SS5. Credits: Esri World Imagery; sources: ESSRI, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community, © Crown copyright 2017. Distributed under the Open Government Licence (OGL). Ordnance Survey (Digimap Licence).

Table 1 Details of the five sub-sites (SSs), which are all currently designated Site of Special Scientific Interest (SSSI), Special Protection Areas (SPA), and Special Areas of Conservation (SAC). FFNNR: Forsinard Flows National Nature Reserve.

3 Materials and methods

3.1 InSAR and time series preparation

InSAR surface motion time series were calculated at a pixel resolution of 80 m × 90 m across the study site. To calculate the surface motion, we used 410 Sentinel-1A and -1B synthetic aperture radar images (descending orbit 125) gathered every 6 to 12 d between 12 March 2015 and 1 July 2019 from the European Space Agency Copernicus Open Access Hub (https://scihub.copernicus.eu, last access: 18 March 2022). Satellite interferometry was applied using these images with the APSIS technique (Sowter et al., 2013). The APSIS technique contains an adapted version of the established SBAS Differential InSAR time series algorithm (Bateson et al., 2015; Cigna and Sowter, 2017). It was designed to improve the density and spatial distribution of survey points to return measurements in vegetated areas, where Differential InSAR processing algorithms habitually struggle due to incoherence (Osmanoğlu et al., 2016; Gong et al., 2016). The APSIS algorithm was implemented using Terra Motion Limited's in-house Punnet software, which covers all aspects of processing from the co-registration of SLC (Single Look Complex) data to the generation of time series (Sowter et al., 2016). Standard interferometry image thresholds were as follows: a maximum horizontal baseline no more than 250 m, a maximum temporal separation of 1 year between image pairs using a coherence threshold of 0.25, and a minimum multi-look point acceptance threshold of 360 to a resolution of 80 m by 90 m. Motion was measured relative to a stable reference point, a building on glacial till, at Wick Airport (58.4533^∘ N, 3.0879^∘ W). Phase unwrapping was implemented using an in-house implementation of the Statistical-cost, Network-flow Algorithm for Phase Unwrapping (SNAPHU; Chen and Zebker, 2001). Two products were produced for each georeferenced pixel location, the multi-annual average velocity (m yr⁻¹) of the time series calculated from the radar line of sight, and the 6–12 d time series of surface motion that detects the seasonal expansion and contraction (or bog breathing) as annual oscillations in relative height (Alshammari, et al., 2018). Each motion time series was then processed as follows.

First, using the R programming environment (R Core Team, 2013), the time series was sub-sampled into equal time intervals of 12 d, to match the longest overpass interval of Sentinel-1 images since Sentinel-1B, which reduces overpass times to 6 d, was not operational until 2016. Outliers were re-estimated using the R “tsclean” function (Box and Cox, 1964), from the R package “Forecast” (Hyndman et al., 2019). Gaps were filled with a linear interpolation using the R “approx” function (Becker et al., 1988) from the R “stats” package (R Core Team, 2020) after “spline” interpolation methods were found to produce contradictory results with adjacent time series across the largest gaps. The “detrend” R function aligned and reset each time series around zero by subtracting the mean.

Second, multichannel singular spectrum analysis (MSSA) using the SSA-MTM toolkit (Ghil et al., 2002; SPECTRA, 2021) was applied to isolate the cyclical, annual seasonal component of the time series from regional climate trends. The MSSA procedure initially calculates covariance after channel reduction with principal component analysis (PCA); then, using moving windows of 2–12 months, long enough to capture annual cycles, we recovered 80% of the signal variance in the first 20 empirical orthogonal functions (EOFs), which included the seasonal cycles in the time series. In the first instance, surface motion time series were reconstructed using EOFs 1–6. This reconstruction captured the seasonal cycles but also included longer-term climate trends, notably 3 wetter years leading to the 2018 European-wide drought (Buras et al., 2020). This climate trend causes merging and shouldering of peaks that compromised the detection of the seasonal cycles, particularly in the west of the study area, where it is wetter. To overcome this difficulty, we used a surface motion time series reconstruction using EOFs 5 and 6 (Fig. 2; Sect 1.1 and Fig. S1 in the Supplement) which extracted only the seasonal cycles. The final MSSA reconstruction provides a signal of relative movement, not absolute surface height, where a rise in the bog surface is an increase in displacement values and a fall is a decrease (Fig. 2).

Figure 2 Examples of surface motion time series and the metric definitions for (a) “soft” wet bog and (b) a “stiff” drier bog, calculated from Sentinel-1 APSIS InSAR time series data between 12 March 2015 and 1 July 2019 against (c) monthly precipitation for northern Scotland (UK Met Office 20-year average, light blue) and Forsinain in the Flow Country (Scottish Environment Protection Agency, Dark blue). The initial mean detrended time series (brown), and multichannel singular spectrum analysis (MSSA) reconstructions (R1) retaining the local climate trend (purple, a combination of empirical orthogonal functions 1–6) and R2 annual seasonal cycles (blue, a combination of empirical orthogonal functions 5 and 6) are shown. The two surface motion metrics used in the analysis are peak timing (t) and amplitude (a) shown for the surface motion year 10 May 2016 to 9 May 2017. A third surface motion metric, multi-annual average velocity, is not illustrated here as it is part of the InSAR data processing (Sect. 3.1). This asynchronous timing of peaks between (a) and (b) forms a bimodal distribution in the peak amplitude timing of the peatland landscape. The drought event of 2018 (Buras et al., 2020) is shaded and can be seen to influence relative surface motion with a local climate trend in the soft, wet bog (a).

Third, the MSSA reconstructions were then analysed for two of three surface motion metrics used to represent the condition of the peatland within each pixel for each motion year using the R “pracma” peak-find function (R Core Team, 2013). Metric one is the date of the annual peak “swelling” in the seasonal cycle (peak timing) of the MSSA reconstruction within 12 months from mid-May (Fig. 2). This has been shown to relate to peatland ecohydrology (Alshammari, et al., 2020; Tampuu et al., 2020). Metric two is the annual maximum amplitude (m) in the surface motion signal (amplitude) measured from the previous minimum of the MSSA reconstruction (Fig. 2). This is an indicator of the elastic response of the peat to changes in water storage (Roulet, 1991; Waddington et al., 2010). Metric three is the multi-annual average velocity (m yr⁻¹) of the peatland surface calculated directly and previously described from the APSIS processing. This is a measure of vertical peatland growth (greater positive value) or subsidence (greater negative value) calculated over a fixed section of the time series (Sowter et al., 2013).

As the 2018 drought caused severe and widespread subsidence, it was found to have subdued the multi-annual average velocity and for this reason, we concentrated our analysis between 10 May 2015 and 9 May 2018, discarding the drought period. Multi-annual average velocity was recalculated accordingly. While the impacts of climate anomalies on the time series were noticeable and interesting, the first step is to gain an understanding of how InSAR data can be used to characterise peatland condition, and we focus on this aspect. We also screened time series with multiple peaks per annum or years where peaks were not discernible, and these pixels were classed as having irregular cycles. Irregular time series made up 8.4 % of the data set and are commonly associated with water courses and damaged bog (including agriculture and some forested pixels). Exclusion of these irregular time series does not affect our conclusions. Additionally, for the first year in the time series (2015 to 2016) many surface motion time series are truncated preventing the accurate calculation of amplitude or peak timing. In that year the mapping is incomplete, so for clarity we show results for a complete motion year from the mid-section of the concentrated analysis (2016 to 2017). To understand the relationships between the three metrics with respect to peatland condition we visualised the metrics in a three-axis plot (Fig. 3).

Figure 3 Characteristics of the motion metrics by sub-site (SS) calculated from MSSA reconstructions of the InSAR-detected annual motion between 10 May 2016 to 9 May 2017. (a) Frequency of peak timing throughout the motion year, where May is the period of least surface motion activity, with a binomial distribution split at 10 November into an early (blue) and late (green) cluster. (b–f) Three-axis plots of the surface motion metrics: (b) SS1 Balavreed, (c) SS2 Cross Lochs, (d) SS3 Knockfin Heights, (e) SS4 Loch Caluim, (f) SS5 Munsary. Axis: x, peak time (Month); y, amplitude (Amp.: m), z, multi-annual average velocity (Vel.: m yr⁻¹). For multi-annual average velocity, greater positive is peatland growth and greater negative is peatland subsidence. Magenta box is for visual reference. Each sub-site demonstrates a specific range in peatland condition according to the plot space it occupies.

3.2 Ecohydrological classification of the sub-sites

The training bed of the sub-sites SS1 to SS5 was divided manually into 130 smaller polygons (hereafter, sub-site polygons). Polygons ranged from (0.3 to 0.6 km²), a practical size (a) for reliable field and map-based validation and (b) to be appropriate for capturing key features of the landscape (e.g. pool systems, forestry or restoration blocks, streams and banks). To find boundaries between the polygons, one of the authors without specialist peatland knowledge searched for distinct contrasts in the landscape structure (e.g. topographic setting, natural drainage, evidence of drainage ditches and where these features would influence hydrology, forestry plantation, restoration, land management and the likely range, consistency or inconsistency in peak timing). This was possible using Google Earth, maps and the InSAR data (Supplement 1.2, Fig. S2). An additional set of 125 polygons was randomly selected across the whole study area using the Esri ArcMap “random point” tool (hereafter, random polygons). The immediate area close to the point were similarly assessed for features in the landscape, to define the polygon boundaries (most between 0.2–0.5 km²). While the sub-sites included the continuum of conditions and features adjacent to each other, the random polygons ensured that there was an improved capture of the ecohydrological states across the whole study area, reducing the likelihood that the sub-sites may have excluded a particular ecohydrological state. Summary statistics of the three surface motion metrics were calculated for each polygon. Average altitude, slope and aspect for all InSAR points in each polygon were also calculated from the Shuttle Radar Topography Mission Digital Elevation Model (Jarvis et al., 2008) to provide measures of topography (Supplement 1.3).

The full set of polygons (sub-sites and random polygons) was then passed to one of the authors with specialist peatland knowledge and based locally for a “blind” (i.e. without prior knowledge of or information about InSAR metrics within the polygons) ground-based ecohydrological classification. For each polygon, the cover of plant functional types (PFTs; Sphagnum, other mosses, shrub, sedges, grasses, rushes, and conifer trees) and the presence of hydrological features (pools, streams, drains, erosion gullies, slope) were recorded using a semi-quantitative scale (0: not present or scarce; 1: present; 2: co-dominant; 3: dominant). Current management (conservation, drainage for agriculture and peat cutting, forestry, restoration by forest-to-bog, restoration by drain blocking, and wind-farm construction) and historical management (e.g. burning, land-use conversion including wind-farm development, and restoration) were also documented for each polygon. The author responsible for classification visited and surveyed all the sub-site polygons and 86 of the random polygons (i.e. 85 % of all polygons) by walking across the area within the mapped polygon. For the random polygons where access was not permitted, shared local knowledge from stakeholders (land managers, project officers on the ground, wardens, and gamekeepers) was used instead of a field visit. In all cases, this was complemented with a combination of existing data, 1:50 000 UK Ordnance Survey maps, NatureScot National Vegetation Classification maps (SNH, 2019), and Google Earth imagery.

Using the semi-quantitative scores, the PFT and hydrology polygon attributes were clustered by similarity using a hierarchical cluster analysis (HCA; Supplement 1.4, Fig. S3). To avoid an overly split hierarchical tree with only one or two members per cluster requiring complex explanation, it was deemed more informative to analyse the PFTs, hydrology, and the topography categories separately from each other. For the PFTs, once the clustering was complete, the average score of the semi-quantitative scale of each PFT in the cluster was ranked. The top three PFTs were used to characterise and name the plant functional group composition (Tables S1–S4 in the Supplement). For data visualisation of the results, clusters were further grouped based on the dominant PFT, resulting in five plant functional groups: Sphagnum, shrub, grass, bare peat (where low or absent vegetation was dominant), and forestry (Table 2). The four polygons that were dominated by rushes (R) had shrub as a co-dominant vegetation, so they were incorporated into the shrub group. While sedges (Sg) were co-dominant in many Sphagnum and shrub polygons, they were not the dominant PFT in any cluster and therefore did not form a separate group.

Table 2 Percentage proportion of clusters derived from hierarchical cluster analysis based on plant functional types (PFTs) represented in the polygons of the five sub-sites and the random polygons. Clusters are defined by the dominant (first) and co-dominant (subsequent) PFTs. PFT notations: Sp – Sphagnum; S – shrub; Sg – sedges; M – moss; G – grasses; R – rushes; F – forest; LoA – low or absent vegetation (brash, bare peat following tree felling or restoration activities, etc.). n: number of polygons. For data visualisation, the sub-site and random clusters were then grouped into five plant functional groups based on the dominant PFTs: Sphagnum, shrub, grass and rushes, bare peat, and forestry. Clusters dominated by rushes (R) were incorporated into the shrub group for data visualisation given their low number and shrub co-dominance. Bare peat and forestry were retained despite low numbers, as their vegetation is associated with specific management intervention. Footnotes describe further differences between clusters (Tables S2 and S4).

^a Additional presence of G. ^b Additional presence of R. ^c Absence of G.

We also categorised topography into equal altitude belts (0–150, 151–300, and 301–450 m) and split slope face direction into four quadrants (north, east, south, and west facing) and ran the HCA. Except for the highest most eroded SS3 site, altitude and aspect did not show any meaningful cluster groups and played no further part in the analysis. The lack of topographic relationships is largely due to the gentle relief of the Flow Country, which has few sheltered slopes and deep valleys. Instead, we used average gradient (degrees) in the polygon and found a natural breakpoint at 1.5^∘ that split the data set equally between flat (<1.5^∘) and slope (>1.5^∘).

3.3 Mapping the state of the peatland system

Within the three-axis plot, we then chose a point with a winter peak timing, a high amplitude, and extreme positive velocity, normally associated with “soft”, wet Sphagnum peat and mapped the whole study area relative to that point. The actual reference point was selected by isolating the points that peaked in the winter and then stepping down through the percentiles of the metrics distributions until the case with the most frequent timing (e.g. February), highest positive velocity (e.g. 0.006 m yr⁻¹), and highest amplitude (e.g. 0.008 m) was identified. For the condition mapping, data for all other points in the three-axis plot were then paired with this reference point and the Euclidian distance between them in three-dimensional (Cartesian) space was calculated. Based on the subsequently observed bimodal distributions (Sect. 4; Fig. 3), if the paired point was on the opposite side of the bimodal distribution to the reference point, the Euclidian distance was mapped as a “V”-shaped path via zero velocity and zero amplitude at the mid-date, 10 November, between the earlier and later timed clusters in the distribution. Prior to calculation, the positions of the outer portions of the three-axis plot were adjusted. This is because if the paired point is before (left of) the earlier peak and after (right of) the later peak the difference between the peak timing and the origin would be overestimated. To mitigate this, these cases were folded inwards along the axis of the peak of their distributions (effectively turning the upturned “W” shape of the bimodal distribution into an “M” shape). Using the natural breaks (Jenks) classification in Esri ArcGIS as a guide, thresholds were used to classify distance from the soft, wet Sphagnum point. These values were then mapped to produce an ecohydrological classification across the whole study site.

To further validate our ecohydrological classification map, we remotely identified and marked the central locations of all the pool systems within the study area (328 in total) using Google Earth images and determined if these markers corresponded to the soft, wet Sphagnum state. Although soft, wet areas in which Sphagnum is dominant do not necessarily contain pool systems, pool systems in this area almost always contain soft, wet peat. Furthermore, in the study area, pool systems provide a spatially distributed, abundant, and easily identifiable sample of this part of the peatland system. They also correspond to the part of peatland systems most unequivocally associated with the “near-natural” ecohydrological condition in this type of upland blanket bog. As the position of each pool marker did not take into account that pool systems often display complex morphology and varying geometries (and sometimes variable condition) related to local hydrological conditions (Goode, 1973; Lindsay, 2018), there was a need to tolerate a level of spatial uncertainty. To capture this, a wider area of 150 m (using the buffer function in Esri ArcMap) was calculated around the marker. This buffer area contained at least three pixels by three pixels of the ecohydrological map. We then calculated the percentage of pixels identified in the ecohydrological classification as soft in the buffer. Pixels classed as irregular were not included in the count (Table S5). We acknowledge that we focus on one particular state of peat, which does not account for the presence of drier, thin, or damaged peat conditions. Indeed, stiffer, thin, and damaged peat cannot readily be associated with a single, well-defined and remotely identifiable feature distributed across the whole study area in the same way. For instance, whilst drains can be observed from Google Earth, their age, maintenance, and the extent of their impact on the peat would require further evidence beyond the scope of this study.

Figure 4 Distribution of mean peak timing date by dominant plant functional clusters for polygons with pools, streams, or other hydrological features (e.g. drains, erosion gullies, peat cutting, or no apparent features) in either a flat (gradient < 1.5^∘) or slope (gradient > 1.5^∘) topographic setting for the random and sub-site polygons. Plant functional group polygons (n random/sub-site): bare – bare peat ($\mathrm{6}/\mathrm{2}$); forestry – conifer plantation ($\mathrm{12}/\mathrm{3}$); grasses and rushes – grass-dominated communities ($\mathrm{13}/\mathrm{14}$), typically Molinia caerulea; shrub – shrub-dominated communities ($\mathrm{59}/\mathrm{77}$), typically either or both of Calluna vulgaris and Erica tetralix; Sphagnum – Sphagnum-dominated communities $\mathrm{35}/\mathrm{34}$. Sedges, rushes, and other mosses are also present, often as co-dominant species in both Sphagnum and shrub communities (see Table 2).

4 Results

From the frequency histograms of maximum peak timing, we discovered a bimodal distribution, showing an early and late peak (Fig. 3a) and defined the motion year to begin at the least active swelling period on 10 May to avoid dividing periods of maximum swelling into consecutive calendar years. The bimodal distribution peaks fall between August to October and December to February, similarly illustrated in the three-axis plots (Fig. 3b–f). For each of the sub-sites, the plots of the three motion metrics show small variations in the shape and position of the data cluster reflecting the diversity of peatland conditions sampled across the landscape.

4.1 Relationship between surface motion and ecohydrology

The HCA revealed ecological groups relating to dominant plant functional types that were comparable between the sub-site and random polygons (Table 2) as well as hydrological groups separating polygons with pool systems and polygons with streams from all other polygons (Fig. S3). When the HCA classifications and topographic information (slope) were compared to the surface motion metrics, we determined the following consistent relationships for the sub-site and random-site polygons.

4.1.1 Timing, hydrology, and topography

Shifts in the peak timing distributions relate to a combination of topography, hydrology, and plant functional group (Fig. 4), and peak timings themselves were consistent within groups between the 3 motion years. Sphagnum-dominated polygons are almost exclusively associated with pools and more so on flat ground with a tendency to have a peak later in the year than most other hydrological and topographical combinations of vegetation or hydrological features. Within the hydrological class pool, polygons with topographic gradients greater than 1.5^∘ (slope) for both random and sub-site polygons have their highest peak timing frequencies earlier in September and October than polygons with topographic gradients less than 1.5^∘ (flat), which tend to peak in November, January, and February. The steeper gradients, at both random and sub-sites, tend to be associated with shrubs and grass-dominated vegetation communities with low frequencies of Sphagnum-dominated polygons. Polygons with pools as the dominant hydrological feature tend to have their highest frequencies from October onwards, later than polygons classed as stream or polygons with other features (drainage ditches, peat cutting, or erosion gullies).

Figure 5 Distribution density of multi-annual velocity for each plant functional group and polygon type plant functional polygons n (random/sub-site): bare – bare peat ($\mathrm{6}/\mathrm{2}$); forestry – conifer plantation ($\mathrm{12}/\mathrm{3}$); grasses and rushes – grass-dominated communities ($\mathrm{13}/\mathrm{14}$); typically Molinia caerulea; shrub – shrub-dominated communities ($\mathrm{59}/\mathrm{77}$), typically either or both of Calluna vulgaris and Erica tetralix; Sphagnum – Sphagnum-dominated communities $\mathrm{35}/\mathrm{34}$. Sedges, rushes, and other mosses are also present, often as co-dominant species in both Sphagnum and shrub communities (see Table 2).

4.1.2 Multi-annual average velocity and dominant plant functional group

Multi-annual average velocities that were most positive (gain of height over time) were almost entirely dominated by Sphagnum (Fig. 5). Polygons with plant functional groups typically associated with natural or man-made drainage (shrub), disturbance (forestry and bare peat), or thin, degraded peat (grasses) consistently displayed negative long-term multi-annual average velocities (loss of mass over time) regardless of topographical setting. Sites in which grasses or forest dominate tend to have a more intermediate multi-annual average velocity than either shrub- or Sphagnum-dominant polygons. Where bare peat is dominant, multi-annual average velocities are lowest (most negative).

4.1.3 Multi-annual average velocity and management

When multi-annual average velocities are compared across different management classes (Fig. 6), the least negative values are observed under conservation management and most negative values are associated with forest-to-bog management, a restoration approach that typically involves compaction from heavy machinery during the removal of conifer stands, followed by drain blocking and surface re-profiling. This restoration class shows a broader distribution in long-term multi-annual average velocity than other management classes, reflecting a variable degree of recovery associated with differing starting conditions, time since initiation (ranging from 0 to >15 years), and techniques used in the intervention.

Figure 6 Distribution density of multi-annual average velocity for different management groups and by polygon type (sub-site and random). Management polygons n (random/sub-site): restoration (forest-to-bog) $\mathrm{12}/\mathrm{13}$, restoration (drain blocking) $\mathrm{0}/\mathrm{6}$, forestry $\mathrm{10}/\mathrm{3}$, drainage $\mathrm{62}/\mathrm{13}$, conservation $\mathrm{37}/\mathrm{95}$.

4.1.4 Time and magnitude of peatland swelling

The factors that influence amplitude can be deduced from relative annual amplitude change and peak timing plots for the three most dominant PFT clusters (Sphagnum, shrub, and grass) across 3 surface motion years (Fig. 7). These graphs all show a positive trend between timing (day of year) and amplitude with a tendency for higher amplitudes to usually occur later in each surface motion year. Steeper slopes tend to have the lower amplitudes that peak earlier in the surface motion year. Shallow to flat slopes tend to have higher amplitudes and peak later in the surface motion year. Year-on-year variation in the range of observed amplitudes is apparent, with a large range in 2015–2016 and smaller ranges in the 2 subsequent years. We attribute this to interannual variation and antecedent conditions in rainfall (e.g. Fig. 2c). A relatively dry 2014–2015 resulted in a strong amplitude response in 2015–2016 with lesser responses in the 2 subsequent wetter years. These differences can be related to the amount of unfilled pore space in the uppermost layer of the peat. In essence, as the peat gets wetter and the pore space fills, there is less capacity in the peat to add more water and the amplitude response diminishes. A more saturated state of the peatland in 2016–2017 and 2017–2018 was also noted in field observations.

Figure 7 Timing of and relative amplitude for 3 consecutive years (2015–2018) with respect to slope gradient (degrees), dominant plant functional type (PFT; grass, shrub, and Sphagnum), and polygon type (sub-site and random).

4.2 Application to large-area condition mapping

The observed relationship between surface motion metrics and ecohydrology is readily interpreted in the context of reported field measurements of peat surface motion (Howie and Hebda, 2018; Morton and Heinemeyer, 2019). Flatter sites under near-natural conditions are poorly drained, wetter, and dominated by Sphagnum spp. In turn, Sphagnum spp. have a considerable capacity for water storage as a direct result of their physiology (Kellner and Halldin, 2002), resulting in peak water storage and seasonal swelling of the surface later in the year. Drier sites with compacted peat have less capacity to store water and reach water holding capacity earlier in the autumn (Price, 2003). Furthermore, the more degraded peat at these sites is less elastic and therefore exhibits a lower-amplitude response to changes in water storage (Holden et al., 2004; Lui and Lennartz, 2019). As the seasonal water balance shifts, drier, better drained sites lose water first followed by the Sphagnum sites, which may continue to swell on account of hysteresis during the first stages of water loss (Howie and Hebda, 2018). Synthesising Sect. 4.2 (Figs. 4–7) the bimodal distribution of peatland surface motion timing within our landscape may be interpreted as reflecting two dominant components of the landscape. Wetter, flatter sites, typically dominated by Sphagnum and sedges, are soft peats; they tend to reach peak surface heights later in the year (December to February window) and have higher amplitudes and more positive velocities. Drier shrub- and grass-dominated sites are “stiff” peats; they tend to reach peak surface heights earlier in the year (August to October window), with slightly (but not exclusively) lower-amplitude oscillations and more negative velocities. Additionally, in the three-axis plot there are points characterised by both low amplitudes, negative velocities, and peak timings that can also fall outside the window of the soft and stiff class. These subtle variations in the metrics identified a third broad ecohydrological class which reflects thin peats, the most degraded and drained grass-dominated sites, or sites under restoration that are in transition to either a soft or stiff state.

Mapping and applying thresholds to the Euclidian distance calculations into these three broad peatland classes relative to the soft peat condition generated the peatland condition map (Fig. 8a). The classification produced a patchwork of conditions in the Flow Country, and the map evidences the widespread occurrence of the stiff peat condition associated with both naturally drier areas on the slopes and wet, soft peat margins but also areas made drier by a land-use history of burning and drainage. The impact of restoration activities, following the felling of forestry on deep peat in the last 25 years, can also be seen: recent forest-to-bog clearance is displayed as “thin/modified” peat, whilst areas well on their way to recovery are showing as either the soft or stiff condition (Fig. 8b and c). In polygons where forest is planted on peat, the signal is much more mixed with a greater proportion of the irregular class. This mixture is a result of the poorer SAR response over trees and the variable conditions encountered within forest stands. For example, in these forestry plantations, fire breaks, gaps between blocks (termed “rides”) and riparian areas that were never planted can still be wet, soft, and Sphagnum-rich in contrast with the planted blocks themselves. Furthermore, plantations may enclose areas of deep peat with pools and may display variable wetness and dryness depending on site and topography.

Figure 8 The ecohydrological classification of “soft”, “thin/modified”, and “stiff” peatland condition, with respect to the location of pool systems, and areas of forest-to-bog restoration in the study area. (a) Classified map based on the Euclidian distance (in three-dimensional Cartesian space) from the position of the soft reference to all other points in the three-axis plot of peak timing, amplitude, and multi-annual average velocity, including the screened irregular time series for the period 10 May 2016 to 9 May 2017. Pool systems that have been correctly classified are shown as a “hit” and otherwise as a “miss”. (b) A detailed view of the classified area highlighted within (a) illustrating the relationship with hydrology, determined from a DEM (potential drainage, PD) and polygons delineating areas for restoration. (c) A true colour satellite image of area showing the restoration status of the polygons, un-felled forest on peat (FoP), peat areas at various stages or different years of forest-to-bog (FtB) restoration, and an area of thin peat (TP). (d) Frequency distribution of Euclidian distance, and the thresholds used to differentiate between the three peat conditions. Credits: the classified area (approx. 930 km²) was delineated using peat soils from the National Soil Map of Scotland (JHI, 2021). Images sourced via Esri ArcMap in 2021. Image source for (a) and (c) Esri World Imagery: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community.

Using the criteria for the remote validation that pool systems should always include a pixel of the soft peat category within the 150 m buffer, 97.9 % ($\mathrm{321}/\mathrm{328}$) of the pool systems were identified as having at least one pixel of that class in the surrounding buffer (Table S5). A detailed inspection of the remaining 2.1 % ($\mathrm{7}/\mathrm{328}$) of pool system markers without a soft classified pixel within the 150 m buffer reveals that these pool systems all showed evidence of localised erosion or drainage causing degradation of their natural hydrology. While not validated in the same way, further inspection of our classification combined with specialist knowledge of the area indicates that the thin/modified class corresponds to areas under restoration, notably areas recently felled for forest-to-bog restoration, areas subject to intensive grazing, and thin peat soils on steeper higher ground and in valley bottoms (e.g. Fig. 8b and c). The abundance of the thin/modified class is striking in the east of the study area, and we note that this corresponds to long-term historical usage of the land for agriculture and associated cutting of peat for fuel (Andersen et al., 2018; Minasny et al., 2019).

Within the area, our method identifies approximately 254 km² (27.3 % of the area) as soft, wet Sphagnum-dominated peat, 481 km² (51.7 %) as stiff shrub-dominated peat, 117 km² (12.8 %) as the thin/modified peat class, and 78 km² (8.4 %) as irregular time series. This classification therefore provides an overall measure of the current state of this blanket bog landscape, to which future regional change, on account of climate change or restoration, may be compared.

5 Discussion

Our most important finding is that surface motion metrics derived from APSIS InSAR time series enable almost continuous spatial and temporal characterisation of peatland condition at large scales. That the SAR data can penetrate cloud cover, measure regular physical displacement of the surface, and capture a known dynamic behaviour associated with peat resilience gives this approach a significant lead over the far more challenging effort to measure peatland condition from optical reflectance data (e.g. Artz et al., 2019). This is compounded by the fact that peatland areas are often obscured by cloud (Minasny et al., 2019). A valuable exercise would be to quantify the similarities and contrasts of motion-mapped peatland condition to optical products, and we anticipate that motion data will bring different and complementary information. This may be advantageous for restoration monitoring, and information on peatland mechanical condition from surface motion may be key to resolving weaknesses in optical studies, for example in carbon accounting (Couwenberg et al., 2011).

The sensitivity and dynamic response of surface motion metrics to changes in the state of the peatland system should make the method ideally suited to monitoring and informing peatland management and restoration. Globally, large areas of northern peatlands degraded by historic drainage, grazing, and forestry are now under or targeted for restoration (Rochefort and Andersen, 2017). Consequently, peatland conservation and restoration are increasingly perceived as critical tools in the fight against global climate change (Leifeld and Menichetti, 2018; Amelung et al., 2020; Günther et al., 2020). Restoration strategies typically involve raising water levels to re-establish wet conditions. The expectation is that this will promote Sphagnum establishment, often a key measure of the success of an intervention (Rochefort and Andersen, 2017; Bellamy et al., 2012; Caporn et al., 2018; González and Rochefort, 2019).

In the case of blanket bog landscapes, our finding of naturally stiff drier shrub and soft, wetter Sphagnum states raises the question as to whether a strategy of increasing Sphagnum cover is always an appropriate restoration target or indeed if it is the only desirable outcome for blanket bogs. In this context, an APSIS InSAR-based assessment of the condition of a whole blanket bog landscape can help guide restoration strategies by first identifying the typical natural states and hydrological structure of that landscape. Second, following intervention, this approach could enable a robust monitoring of restoration trajectories and outcomes (Marshall et al., 2021). Our study also suggests that when monitoring restoration trajectories over time, the impacts of interannual variabilities such as precipitation on the metrics would likely need to be accounted for.

In natural landscapes, these peatland states are a consequence of landscape evolution in which the vertical accumulation of peat must be counterbalanced on an appropriate spatial and temporal scale by erosion (Large et al., 2021). Drier states correspond to areas of net carbon loss due to natural drainage, incision, and erosion along peatland margins, and wetter states correspond to peatland interiors, areas with low gradient that tend towards carbon accumulation. In this context, to restore a site that is naturally stiff and dry to the soft, wet state would risk instability, while the opposite would fail to optimise carbon storage. A more suitable and sustainable ambition is to accept that restored blanket bog sites may follow different trajectories towards naturally Sphagnum or shrub states (or soft and stiff) and that these target end states will be constrained by the hydrological landscape setting, as conceptualised by Winter (1988). Our approach provides evidence for these natural states co-existing within the study areas and evidence to guide and monitor appropriate restoration trajectories within this system. Recognising and preserving this mosaic is critical in maintaining large- and small-scale peatland landscape stability and carbon balances, particularly as long-term models suggest that the natural drying out of peatland is accelerating due to drainage (Harris et al., 2020; Leifeld et al., 2019) and climate change (Gallego-Sala and Prentice, 2013).

The approach outlined here should be readily transferable to alternative peatland settings within different parts of the global peatland climate space. Using surface motion metrics identified from the InSAR time series of peatland motion, a surface deformation space for a given peatland system can be defined. The position of ecohydrological characteristics within this space can then be deployed to quantify the state of the peatland system and map changes with respect to climate change and management intervention. This capacity to customise the approach is valuable as it provides the means to measure peatland condition at a global scale. If realised, this would enhance our understanding of the large-scale functionality of peatland landscapes and provide the robust evidence base required for sustainable peatland management.