Long-term trajectory and temporal dynamics of tropical peat subsidence in relation to plantation management and climate

long-term trajectory of subsidence, are necessary to determine the future economic and environmental sustainability of managed peatland landscapes. While many peatlands in temperate regions such as Europe have been drained for centuries, those of Southeast Asia have mostly been drained for agriculture and forestry practices within the last 30 years. These areas are subsiding rapidly, but few long-term subsidence records exist, and it is unclear whether currently high subsidence rates will be maintained in future. Furthermore, large-scale climate fluctuations associated with the El Ni ˜ no Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD) exert a strong year-to-year influence on rainfall rates, affecting water table depth dynamics in both managed and intact peatlands. In this study, we analysed data collected over more than a decade of subsidence measurements from over 400 plantation and forest plots in Sumatra, Indonesia, including a unique set of 62 sites that have been continuously monitored since 2007. We show that spatial and temporal variations in subsidence rates are primarily determined by water table depth. We also find evidence of declining subsidence rates as a function of time since initial drainage, consistent with previous instrumental records from high-latitude peatlands and recent satellite data from tropical peatlands. Subsidence rates over the study period were strongly affected by the large ENSO/IOD-linked drought event in 2015 – 16, which caused an acceleration of subsidence across all sites. In plantation areas, we estimate that this climate perturbation caused around 14 % of subsidence observed over a twelve year period. At interior forest sites this rose to 32 %, and we found little evidence of ecosystem recovery to the end of 2018. This raises the possibility that repeated extreme droughts in the region could lead to long-term degradation of peat swamp forest ecosystems.


Introduction
Subsidence of drained peatlands has been a recognised challenge for agricultural and environmental sustainability since at least the mid-20th century (Thompson, 1957).Many areas that have been subject to longterm drainage, including large parts of the Netherlands, Northern Germany, Eastern England, the San Joaquin Delta in California, and Southern Florida, have subsided by several metres, to the extent that large areas of agricultural land are now below sea-level, and require protection from flooding by embankments and pumped drainage.Contemporary subsidence rates in these high-latitude peatlands are now typically 1 to 3 cm yr − 1 (Evans et al., 2019;Hoyt et al., 2020).This subsidence is the result of a combination of ongoing soil compaction following the lowering of water levels, and continuous oxidation of aerated peat above the water table.Peat oxidation leads to CO 2 emission from peatlands drained for agriculture and forestry, making them significant contributors to global greenhouse gas (GHG) emissions in many countries (Leifeld et al., 2020).
In Southeast Asia, most peatland conversion has occurred more recently, with around 60 % of the tropical peat swamp forests of Malaysia and Indonesia having been converted since 1990 (Wijedasa et al., 2018).A total peatland area of 78,500 km 2 has been developed for agriculture and silviculture, of which 43 % is under cultivation by smallholders, 39 % in industrial oil palm plantations, and 11 % in pulpwood (Acacia crassicarpa) plantations (Wijedasa et al., 2018).All of these activities involve some degree of peat drainage, and therefore lead to subsidence.In the period immediately following land conversion, rapid 'primary' subsidence, which results from consolidation of the peat after de-watering, can be in the order of a metre within the first two years (Andriesse, 1988;DID and LAWOO, 1996;Hooijer et al, 2012).Subsequent 'secondary' subsidence rates in drained tropical peatlands are typically in the range 2-6 cm yr − 1 (Evans et al., 2019;Hoyt et al., 2020).Higher rates compared to high-latitude peatlands are partly attributable to the higher temperatures experienced by tropical peatlands, which lead to high rates of oxidation when peatlands are exposed to aerobic decomposition, and may also be related to the relatively short time since drainage.As a result, CO 2 emissions from drained tropical peatlands are estimated to be approximately double those from highlatitude peatlands exposed to a similar degree of drainage (IPCC, 2014;Evans et al., 2021).Since many drained Southeast Asian peatlands are in low-lying coastal areas, it has also been suggested that ongoing subsidence could lead to drainability issues in plantations, and to future flood risk (Sumarga et al., 2016;Saputra, 2019;Lupascu et al., 2020).
There is however an ongoing debate as to whether high observed subsidence rates in recently drained tropical peatlands will be maintained over the longer-term.In Northern peatlands, for which several very long-term (>100 year) records exist, subsidence rates typically follow an exponential decay curve under conditions of constant drainage (Hutchinson et al., 1980;Stephens et al., 1984;Deverel and Leighton, 2010).This is attributable to a combination of increasing mineral content and decreasing porosity in the residual peat, which reduces compaction, oxygen transport and thus the amount of organic matter exposed to aerobic decomposition, and to the decreasing volume of peat above the water table, which again reduces the organic matter pool exposed to decomposition.Rapid subsidence can however be re-initiated if drainage levels are increased in order to maintain agricultural productivity (Hutchinson et al., 1980;Stephens et al., 1984).
In tropical peatlands, the lack of equivalent records precludes a truly long-term assessment of subsidence trajectories.Some studies have predicted that tropical peatlands will follow a similar trajectory to highlatitude systems, with subsidence rates gradually slowing over time (e.g.Muruyama and Bakar, 1996;Wösten et al., 1997;Othman et al., 2011), a finding that appears to be supported by recent analysis of satellite data (Hoyt et al., 2020;Umarhadi et al., 2022), but which has yet to be confirmed by long-term ground observations.Others have argued that subsidence in managed plantations will proceed at a more or less linear rate following the initial period of rapid subsidence (Hooijer et al., 2012;Couwenberg and Hooijer, 2013).This argument is based, in part, on the very low ash content (usually less than 2 %) of most tropical peats, which means that there is little or no 'dilution' of the aerated organic matter pool by accumulating residual mineral material.Couwenberg and Hooijer (2013) also point out that in highly managed plantation landscapes, water levels in drainage canals will effectively 'follow the peat surface down' in order to maintain constant relative drainage depths and maintain crop yields, resulting in a constant depth of aerated peat, and near-constant subsidence until the water table intersects the underlying mineral substrate.On the other hand, recent regulatory developments, notably in Indonesia, require that water levels are held at higher levels within plantations, which has the potential to slow rates of subsidence (Evans et al., 2019).
Improved understanding of the long-term trajectory of subsidence in cultivated tropical peatlands is of high importance from both an economic and an environmental perspective.If subsidence rates are indeed linear, then the onset of drainability and flood risk issues may be relatively rapid.If on the other hand subsidence rates decrease over time, these issues may be delayed, potentially by decades.The linearity or otherwise of subsidence also provides an indirect indication of wheher or not current levels of observed CO 2 emissions are likely to remain constant in future.However, because current projections of future subsidence are reliant on short time series, model assumptions are hard to validate against available data.Furthermore, the period over which recent subsidence data have been collected in Southeast Asia has been affected by large climatic perturbations, most recently the strong combined El Niño Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD) event of 2015-16 (hereafter referred to for simplicity as the 2015-16 El Niño) which led to a severe dry season in the region.As a result, it is difficult to disentangle the relative importance of drainagerelated management from that of variations in climate.It is also difficult to determine the extent to which accelerated subsidence during drought periods may be considered reversible (i.e. a broadly natural fluctuation superimposed on the drainage-related trend, for example due to short-term shrinkage and swelling of the peat) or irreversible (i.e. an acceleration of existing drainage impacts, for example as a result of additional oxidation).These potential interactions may also differ between plantations with managed water levels, and forests without direct management impacts.
In this study, we analyse temporal variations and the long-term trajectory of subsidence in a uniquely long-term (12 year), large-scale and methodologically consistent dataset of subsidence measurements obtained from one of the largest tropical peatland regions in Indonesia, Riau province in Sumatra.The analysis follows a previous study (Evans et al., 2019) which quantified average rates of subsidence in industrial Acacia plantations and adjacent conservation-managed native forests, and evaluated the factors influencing spatial variations in subsidence rate.This analysis, based on 312 subsidence monitoring points for which at least three years of data were available, showed an average subsidence rate of 4.3 cm yr − 1 under Acacia, with similar rates of subsidence in forest edge locations, declining to a lower (but still substantial) average subsidence of 2-3 cm yr − 1 at forest interior sites (i.e.those > 300 m from the plantation edge).The reasons for this apparent subsidence at forest interior sites was unclear, with possible explanations including an extended impact of plantation drainage on forest hydrology, a more widespread degradation of peat domes within the region (for example due to large-scale climatic changes or to historic disturbances such as logging activity), or the impact of the 2015-16 El Niño.In the current study, we first analyse drivers of spatiotemporal variation in the full dataset of 434 subsidence poles, and then analyse in detail a subset of 62 measurement sites that have been in continuous operation since monitoring began in 2007.Our aims were to: i) examine the spatiotemporal coherence of subsidence variations at different sites; ii) evaluate the long-term trajectory of subsidence in both plantation and native forest areas; iii) disaggregate the effects of drainage and ENSO/IODrelated climate perturbation on overall rates of subsidence within plantation, forest edge and forest interior areas.Based on the results, we consider the implications of different subsidence trajectories for the management and conservation of tropical peatlands.

Study area and measurement sites
The coastal lowlands of Eastern Sumatra hold one of the world's largest contiguous areas of tropical peatland (Vernimmen et al., 2020).Peat largely formed during the last 8000 years, and extends over C.D. Evans et al. approximately 70,000 km 2 (Dommain et al., 2014).Peat thickness exceeds 3 m over almost all of the area, and exceeds 10 m in some areas, while surface elevations rarely exceed 10 m above sea level.Mean annual rainfall at Pekanbaru Airport, Riau Province, is 2900 mm.There are two wet seasons, during which mean monthly rainfall exceeds 300 mm, in November-December and March-April.There is a short dry season in February, and a longer dry season from June to August, during which mean monthly rainfall falls below 200 mm.Air temperatures are high and stable, with monthly means ranging from 29 to 32 • C (Badan Meteorologi, Klimatologi dan Geofisika, 1994-2017 data).
We analysed subsidence data collected by Asia Pacific Resources International ltd (APRIL) in Riau, Sumatra, from Acacia plantations and adjacent native forest areas.The study area encompasses Measurements includes the Kampar Peninsular, a 6800 km 2 peat dome with a large intact central forest area, surrounded by plantations, and other large peatland areas including the islands of Pulau Padang and Pulau Rupat to the north, and areas adjacent to the Kerumutan National Park to the south (Fig. 1).The Acacia plantations were established from 1992 onwards, and are managed on a five-year rotation from planting to harvest.Peat in plantation areas is mainly fibric (poorly decomposed with abundant plant remains) but becoming hemic (fewer visible plant remains) towards the surface (Hooijer et al., 2012).In undeveloped forest areas, the majority of the peat is fibric.To support plantation growth, water table depths (WTDs) are managed within topographically-defined water management zones, comprising large canals either side of a pair of plantation compartments, and smaller field drains within each compartment, with water levels within each zone controlled by outlet sluices (see Evans et al., 2019 for further details).Until recently, Acacia plantations were managed for a target WTD of 70 cm.However the El Niño fires of 2015 led the Indonesian government to establish a Peat Restoration Agency (Badan Restorasi Gambut, BRG) with a remit to reduce fire incidence and restore 2 million ha of degraded peatland by 2020, and to introduce new regulations (most recently SK.22/PPKL/PKG/PKL.0/7/2017) requiring that water tables be maintained within 40 cm of the peat surface at the centre of each plantation block for at least half of the year, and within 100 cm of the Fig. 1.APRIL concession areas on peat and subsidence monitoring locations, Riau, Indonesia.Blue circles indicate subsidence poles that have been monitored since 2007 (used for time series analysis), red circles indicate sites that were established more recently (also included in the multivariate analysis).Land cover is shown for APRIL concession areas only.
C.D. Evans et al. surface at all times.No active management of water levels occurs within the forests, but there is evidence that peripheral forest areas (within 300 m from the last canal in the plantation) are affected by water management in the adjacent plantations (Evans et al., 2019).For interior forest sites (>300 m from the nearest canal) the same analysis showed no evidence of a direct drainage impact, although water levels appear generally lower than would be expected in a truly pristine peat swamp forest, for reasons that remain uncertain (Evans et al., 2019;Deshmukh et al., 2021).

Subsidence and water table depth measurements
Since 2007, APRIL has monitored an expanding number of subsidence poles throughout its plantation and conservation forest concessions in Riau (Fig. 1).For this study we analysed data from a total of 434 poles -319 in Acacia plantations, 115 in conservation (native) forest.These sites extend across much of APRIL's production and conservation concession areas.We then undertook a more detailed temporal analysis of data from a subset of 'long-term' subsidence poles that were installed at the start of the monitoring programme in 2007, and which have continued (with few or no missing quarterly data points) until December 2018.This subset comprised 62 sites, of which 38 were in Acacia plantation, and 24 in conservation forest.The sites were located along nine monitoring transects, of which seven are in the western part of the Kampar Peninsular, and two in the smaller Mandau plantation area to the northwest.Four of the transects provided data from plantation sites only, three from forest sites only, and two from both plantation and forest sites.The forest sites were located in forest fragments within the plantation area, within buffer zones (areas of native forest of approximately 300 m width on the boundary ofbut withinthe plantation concessions), and in adjacent conservation-managed concessions.The most remote long-term site was 1.4 km from the nearest plantation canal.
At each location, peat subsidence is measured using a hollow, perforated 5 cm diameter PVC pipe, inserted vertically into the peat and anchored into underlying mineral subsoil (Fig. 2).Ground elevation change is recorded quarterly by measuring distance from the top of a pole (the datum point) to the ground surface.Concurrently, WTD is measured by recording the distance from the pole top to the water table, and subtracting the distance from the pole top to the ground surface.While subsidence represents a cumulative measure of peat elevation change since the previous measurement, the WTD measurement is simply a 'snapshot' of conditions at the time of sampling.Further details of the site network and measurements are provided in Evans et al. (2019).

Supporting meteorological data
Long-term records of rainfall were obtained from a network of 22 weather stations operated by APRIL across their concessions throughout the monitoring period.To provide a single, spatially-averaged index of rainfall across the study area we took a mean of the total measured rainfall across all sites during each quarterly period.We also obtained data for the ENSO Index (https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php) and the IOD Index (https://psl.noaa.gov/gcos_wgsp/Timeseries/DMI/).These indices describe east-west sea surface temperature gradients in the equatorial Pacific and Indian Oceans, respectively.The ENSO represents one of the largest climate cycles in the Earth System, affecting weather conditions in many parts of the world.In Indonesia, positive ENSO Index values (i.e.El Niño events) are associated with higher temperatures and extended dry seasons, contributing to ecological degradation of peatland regions through increased incidence and severity of fires (van der Werf et al., 2008;Nechita-Banda et al., 2018).The IOD is a linked phenomenon occurring within the Indian Ocean, with the potential to exacerbate climate perturbations in Sumatra (Saji and Yamagata, 2003).All of the meteorological and climate index data were aggregated up to quarterly means or totals, to correspond to the frequency of subsidence measurements, and to annual (calendar year) values.

Full-dataset multivariate analysis
To evaluate the combined effect of a wide range of locally measured environmental variables on subsidence rates, we implemented a multivariate regression analysis on the full dataset of 434 subsidence poles.The analysis included variables which vary both spatially and temporally (e.g.WTD), temporally only (a single rainfall time series was applied for all sites), and spatially only, including both continuous variables (e.g.distance from canal) and categorical variables (e.g.tree species).Variables used are shown in Table 1.TIME refers to time (in quarters) since the start of monitoring.DCANAL represents the distance (m) of the measurement point from the nearest plantation canal, while DFOREST represents the distance (m) to the nearest forest edge (plantation sites only).PLANT and FOREST are dummy variables (1/0) describing site type.Within the plantation areas, sites were assigned a value of 1 for the dummy variable ACACIA if they were located in Acacia crassicarpa stands, or OTHER if a different plantation species was present (e.g.Melaleuca sp.).Because all sites are either plantation or forest, we excluded FOREST from the multivariate analysis, and similarly for tree species within plantations we excluded the variable OTHER.Peat thickness was not included in the analysis, because previous work on the same dataset suggests that it has little influence on subsidence rates (Evans et al., 2019), and because around 20 % of the dataset lacked peat thickness data.We also lacked sufficient information on peat properties at all subsidence poles to enable detailed peat classification (e.g.Veloo et al., 2014) so were unable to include this explicitly in the analysis.However as noted above, peat is consistently fibric in forest areas, whereas it grades from fibric at depth to hemic near the surface in the plantations, so broad effects of peat type could be captured by the PLANT summary variable.
Prior to analysis, we removed data outliers using an objective procedure to avoid anomalous measurements exerting undue influence on the results.According to Baltagi (2011) such influential observations can be attributable to: (i) improperly recorded data; (ii) observational errors in the data; (iii) misspecification of the statistical model (i.e.wrongly included or excluded explanatory variables); or (iv) legitimate outlying data points.The subsidence and WTD data were collected for operational rather than research purposes, and following consultation with APRIL field staff and initial assessment of the data, we considered that outliers were most likely to be genuine errors (e.g.due to recorder error or site disturbance).We employed the blocked adaptive computationally efficient outlier nominator (bacon) procedure in Stata to detect and then exclude these outlying observations from our dataset.The original algorithm of the bacon procedure was proposed by Billor et al. (2000), and enhanced by Weber (2010).The procedure implements the Mahalanobis Distance approach, which has been widely recognized in statistical analysis, and used in a previous analysis of this dataset by Evans et al. (2019) to detect and exclude outliers for calculating average subsidence rates.Outlier exclusion reduced the total number of observations used in the analysis, as shown in Table 2. Outlier exclusion reduced the range and standard deviation of SUBS, but had relatively little impact on mean SUBS, or on the other variables.
We also checked for correlations between explanatory variables, to avoid problems of collinearity.The strongest correlation was between WTD and RAIN (r = 0.37), well below the threshold of 0.75 often used as a cut-off (Hair et al., 2010), so we retained all variables in the analysis.Because some explanatory variables were only applicable to plantation or conservation forest areas, we applied three regression equations: Plantation and Forest (1); Plantation only (2); and Forest only (3).Each equation was regressed using Pooled Least Square (PLS) specification with robust standard error.
All sites: Plantation sites: Forest sites: In addition to the variables already specified in Table 1, α is the intercept, the β terms are coefficients, and ε is the residual.The subscript i indicates pole, and t indicates time in quarters.Subscript × indicates lags.For SUBS, the lags are from 1 to 4 (i.e. from the preceding quarter up to one year, to account for autocorrelation).For WTD and RAIN, the lags are from 0 to 4 (i.e. from the time of subsidence measurement up to one year before).

Analysis of long time series
Whereas the multivariate analysis above was used to explore the overall drivers of spatio-temporal variation across the entire dataset of > 7500 individual observations, the time-series analysis aimed to provide more detailed understanding of the temporal trajectory of subsidence as a function of plantation management and climatic variation, using data from sites with complete 12 year time series.Based on information provided by APRIL, we assumed that there was no long-term change in drainage intensity at individual monitoring sites over the monitoring period.This reflects the management of the plantations (prior to new Indonesian Government regulations) to achieve a target mean WTD of 70 cm, via adjustment of sluice levels at the outflow of each water management zone.This assumption does not preclude shortterm fluctuations in WTD in response to periods of high or low rainfall, or spatial variations in the mean WTD at different sites due to local variations in topography or proximity to a canal.Adjustments in drainage depth made to meet the new Indonesian water table regulations, which came into force in 2018, were assumed to have occurred too late in the monitoring period to influence results.
Based on the large-dataset analysis of Evans et al. (2019), which showed no clear spatial controls on variation in subsidence within the plantation dataset, and declining subsidence with distance from the nearest canal within the forest dataset, we divided our time series dataset into three subsets, namely: 1) all plantation sites; 2) 'forest edge' sites, located<450 m from the nearest canal; 3) 'forest interior' sites located>450 m from the nearest canal.The choice of cut-off value within the forest sites was partly informed by the previous analysis, which suggested stronger plantation impacts within 300 m of the plantation boundary, and partly pragmatic as the long-term forest subsidence data split naturally into two similarly sized subsets (13 sites < 404 m to the nearest canal, 11 sites > 504 m to the nearest canal).
Data from each land cover category were then analysed according to the following procedure.Firstly, we subtracted the initial pole height from each individual time series in order to derive a time series of elevation change relative to an initial value (in December 2007) of zero.Sites within each subset were then aggregated by taking the median, 10th and 90th percentile relative elevation value for each quarterly measurement interval.This analysis provides an indication of the average trajectory of change within each category, as well as the between-site variability in absolute rates of subsidence.However, given the small size of each dataset this approach is sensitive to missing values (over the full dataset 5.4 % of quarterly measurements were missing), and differences in absolute rates of subsidence will to some extent mask underlying coherence in temporal behaviour.Therefore, we also transformed each individual time series into consistent dimensionless units ('Z-scores') by subtracting the mean of the dataset and dividing by the standard deviation (Equation ( 4)): where x is the observed value, μ is the mean of the data for that location, and σ is the standard deviation.This transformation gives all time series a mean of zero and a standard deviation of 1.For each quarterly mea- surement interval, we then calculated the median, 10th and 90th percentile Z-score for the full set of sites in each data subset.This approach, which has been widely applied to the analysis of multi-site water quality monitoring datasets (e.g.Watmough et al., 2004;Davies et al., 2005;Evans et al., 2010;Oulehle et al., 2013) provides an effective measure of the degree of spatio-temporal coherence within the dataset.High-coherence datasets will have narrow 10th-90th percentile Z-score ranges, whereas low-coherence datasets will have wide ranges.
The approach also provides an indication of the relative importance of short-term variability relative to the long-term trend; datasets dominated by long-term trends will exhibit large changes in median Z-score over the monitoring period, whereas datasets dominated by short-term fluctuations will have large inter-year variations in median Z-score compared to the full-dataset trend.For a more detailed description of this approach, see Evans et al. (2010).
For the quarterly WTD measurements, which were made concurrently with the subsidence measurements, we applied an identical standardisation procedure.For the rainfall, ENSO and IOD index data, we calculated annual and quarterly means.Relationships between hydrometeorological variables and subsidence were analysed by linear regression with the residuals of subsidence trends, described below.

Attributing subsidence rates
We fitted three alternative best fit lines to the median Z-score subsidence rates for plantation (Z pl ), forest edge (Z fe ) and forest interior (Z fi ) (Equation 5).First, we fitted a simple linear regression to the full dataset (Z pl all , Z fe all , Z fi all ).Since the full measurement period includes the 2015-16 El Niño, this model incorporates both land-use and climaterelated contributions to subsidence.To exclude the potential influence of the 2015 event, we also fitted a linear regression to median Z-score data collected before the start of 2015 (Z pl pre , Z fe pre , Z fi pre ).This approach provided an estimate of linear subsidence in the absence of large climatic disturbances, and could be considered an estimate of the average rate of subsidence that is attributable to land-use activities.
where Z i is the mean Z-score in dataset i = [Z pl all , Z fe all , Z fi all, Z p pre , Z fe pre , Z fi pre ], a i and b i are regression coefficients, t is time elapsed from the beginning of observations (January 1st 2007, t in years) and ε is the residual term.Finally, to allow for possible non-linear rates of subsidence, we also fitted an exponential curve (first order decay model) to the same pre-2015 dataset, of the form: (-kt)  (5) where Z i represents the Z score of median subsidence for each landcover class i as above, t is time in years since January 1st 2007, and m and k are fitted constants (m represents the starting value and k the decay rate in yr − 1 ).An arbitrary value of 10 was added to Z i to ensure that all values were positive during both the observation period and the projection period.This approach allows for a gradual slowing of subsidence rates over time, and also excludes the effects of large climate perturbations.

Hydrometeorological data
The 2015-16 El Niño event is clearly evident as a prolonged and strongly positive series of both climate indices, with peak IOD index values following the ENSO index peak (Fig. 3).There was a smaller ENSO index peak in 2009, with negative values for most of the remainder of the study period.The IOD index showed a similar small positive phase in 2009-2010, followed by a more extended positive phase rising from the end of 2012 to a peak in early 2016.On both occasions the IOD index remained positive after the ENSO index had returned to neutral or negative values.
Rainfall data show seasonally alternating wet and dry periods.The third quarter of 2015, during the El Niño event, was the driest threemonth period in the record, and the following wet season was the weakest recorded.The weaker 2009 El Niño event led to two quarters of below-average rainfall, but conditions were not exceptional compared to other time periods.Water table depths were also only slightly affected by the 2009 El Niño, and remained predominantly above the wholeperiod mean (i.e.positive median Z score) until the end of 2012, when both the ENSO and IOD indices were mostly in a negative phase.From 2013 onwards, when the IOD index was consistently positive, water levels were predominantly below the long-term mean, while the strong 2015-16 El Niño led to extremely low water levels over a sustained period through 2015-16.During 2017-2018, when the ENSO and IOD indices returned to mainly negative values, WTDs were similar to  Over the full record (based on the subset of 62 sites with continuous records) median quarterly water table Z scores showed a weak positive correlation with quarterly rainfall totals (R 2 = 0.16, p = 0.009), a weak negative correlation with the IOD index (R 2 = 0.13, p = 0.014), and a stronger negative correlation with the ENSO index (R 2 = 0.28, p < 0.001).Comparing water table Z scores for the three land cover categories shows a high degree of temporal coherence (Fig. 4).Median water table Z scores were strongly correlated between the three landcover categories (plantation vs forest edge R 2 = 0.60, plantation vs forest interior R 2 = 0.46, forest edge vs forest interior R 2 = 0.71, all p < 0.001).

Full-dataset regression analysis
Table 3 shows results of the multivariate regression analysis for the full dataset, the plantation and forest subsets.For all three regression equations, WTD, RAIN and TIME and their lags were in general highly significant (p < 0.01) in explaining spatiotemporal variations in subsidence, with mostly negative coefficients indicating that statistically quarterly subsidence was lower compared to subsidence in previous periods, showing general pattern of lower subsidence rates over time.Meanwhile, in general the results also indicate that subsidence was lower when water tables were shallower, rainfall was higher, and more time had elapsed since the start of monitoring.The fitted coefficient for lagged SUBS was highest for the forest subset.Fitted coefficients for WTD and RAIN were also highest for the forest subset.The magnitude of coefficients appears physically plausible; for example, raising water tables by 10 cm in plantation areas would reduce annual subsidence by 0.86 cm yr − 1 .The effects of time were strongest in the plantation dataset, with subsidence reducing on average by 0.084 cm yr − 1 , compared to 0.061 cm yr − 1 in the forest dataset.The variable PLANT was also strongly significant and negative for the full dataset, suggesting that subsidence rates are on average 0.80 cm yr − 1 lower in plantation versus forest sites if all other conditions are equal.
For the plantation-only dataset, the spatial variable DCANAL was also significant, indicating that subsidence is reduced by 0.36 cm yr − 1 for every 100 m from the nearest canal.Other variables describing distance to the nearest forest and tree species (Acacia versus other species) were not significant.For the forest-only dataset, the significant predictors were antecedent SUBS, WTD, RAIN and TIME (i.e.distance from the nearest drainage canal did not explain any additional variance beyond the covarying effect of WTD).The coefficient for WTD was similar to the plantation-only and full-dataset analyses, but the coefficient for RAIN was larger (-0.0023 for the forest-only versus − 0.0006 for the plantation-only dataset), suggesting that subsidence rates in the forest are more responsive to rainfall variation.The R 2 for the forest dataset was the highest overall (0.192 versus 0.140 and 0.131 for the full-and plantation-only dataset respectively).

Subsidence time series
The Acacia plantation data show fairly steady subsidence over most of the monitoring period (Fig. 5a), but with a 'step change' quarterly subsidence of almost 10 cm between September and December 2015, which was followed by slower average subsidence during 2017-18.The 10th to 90th percentile range increased steadily over the monitoring period, reflecting different rates of subsidence at different sites (Fig. 5).Overall, plantation sites recorded a median subsidence of 51 cm (10th-90th percentile range 32 to 82 cm) over the 12 year monitoring period.
Forest sites within 450 m of a plantation canal showed a similar pattern to the plantation sites, with steady subsidence until 2014, high subsidence during 2015, and limited subsidence thereafter (Fig. 5b).Median subsidence in this dataset over the 12 year monitoring period was 35 cm, with a 10th-90th percentile range of 22 to 64 cm.Forest sites>450 m from the nearest plantation canal (Fig. 5c) showed little or no subsidence until 2010, followed by approximately linear subsidence until the end of 2014, a sharp decrease in 2015-16, and little further change during 2017-18.Median total subsidence for this dataset was 27 cm, with a 10th-90th percentile range of 17-35 cm.
Standardisation of the data greatly reduces the range of apparent between-site variability (Fig. 6), indicating a high degree of temporal coherence in underlying subsidence dynamics within each land-cover category, regardless of the absolute rate of change.This is most strongly evident for the plantation sites, as reflected in the narrow 10th-90th percentile ranges (mean 0.34 standardised units) compared to the change over the monitoring period (around − 3.16 standardised units).For the forest edge sites, the mean 10th-90th percentile range is slightly wider at 0.54, and the overall change slightly smaller at − 3.02.For forest interior sites the equivalent figures are 0.52 and − 2.53.The forest

Table 3
Fitted coefficient values for full-dataset regression analysis of drivers of spatiotemporal drivers of subsidence (in cm quarter -1 ), see Equations 1-2 Table 3. Fitted coefficient values for full-dataset regression analysis of drivers of spatiotemporal drivers of subsidence (in cm quarter -1 ).The dependent variable is SUBS (t), independent variables are antecedent quarterly subsidence from the preceding quarter (t-1) to one year (t-4); water table depth and rainfall at time t and for each antecedent quarter.sites thus show higher internal variability relative to the long-term trend.Nevertheless, strongly coherent behaviour is evident in all three categories, and all show a clear acceleration of subsidence during the 2015-16 El Niño.This is most evident at the forest interior sites, where around half of all observed subsidence occurred during this one-year period (Fig. 6c).

Empirical models of long-term subsidence
The three alternative regression model fits to the median Z score values for each land-cover class are shown in Table 4 and Fig. 6.All three regression models were highly significant for all three land-cover categories, although they were consistently weakest for the Forest Interior sites.The full-period linear regression model gave subsidence rates of − 5.17 cm yr-1 in the Acacia plantation, − 3.76 cm yr-1 in the forest edge sites, and − 2.86 cm yr-1 in the forest interior sites.The linear model fitted to the pre-2015 data gave subsidence rates of − 4.46, − 3.49 and − 1.86 cm yr-1 respectively.The effect of including or excluding the 2015-16 El Niño event in the linear model was smallest for the forest edge sites, and greatest for the forest interior sites (compare black and red lines in Figure 6).
The pre-2015 exponential model gave marginally stronger correlations than the linear model fitted to the same period for the plantation and forest edge datasets, but a slightly weaker correlation for the forest interior dataset, consistent with the inclusion of time since drainage as a predictor in the multivariate analysis (Table 3).The model indicated a reduction on the subsidence rates over 12 years monitoring period (excluding ENSO/IOD impacts) from 4.89 to − 3.74 cm yr-1 in the Acacia plantation, − 3.85 to − 2.91 cm yr-1 in the forest edge sites, and − 2.01 to − 1.67 cm yr-1 in the forest interior sites.While differences in model performance are slight (Table 4) and there is very limited divergence between linear and non-linear models fitted to the pre-2015 data (Fig. 7, comparing green and red lines), the exponential model is the most consistent with the results of the multivariate analysis in showing a negative effect of time on subsidence rates.

Influence of the 2015-16 El Niño
To evaluate the influence of the El Niño event on post-2015 subsidence, we calculated residual Z scores for subsidence in each land cover category after removing the linear subsidence trend for the pre-2015 period (red line in Fig. 6) and compared these to median water table Z scores in the same category (Fig. 7).In all cases, we found that the El Niño led to a 'step change' drop in peat surface elevation between September and December 2015.Furthermore, for all categories we observed a significant linear correlation between residual subsidence Z score and water table Z score during both the pre-and post-El Niño periods (Table 3), but in all cases this relationship was offset after the El Niño event (compare grey and black circles in Fig. 7).At the plantation sites, the gradient of the peat elevation versus water table relationship was similar before and after the El Niño, while for the forest edge sites it became steeper.Conversely, in the forest interior sites the relationship between peat elevation and water table (which was very strong before the El Niño event) became markedly weaker afterwards.
Finally, we plotted quarterly residuals of the 2007-2014 linear regression model against time in Fig. 8.The data show little deviation from the linear trend during 2007-2014, but a clear and consistent negative deviation for all land-cover categories from 2015 onwards.For the plantation and forest edge sites, residual Z scores were close to zero in the most recent measurements, suggesting that peat elevations have returned to the levels that would have occurred in the absence of the El Niño.This was not the case in the forest interior sites (>450 m from the nearest canal), where residuals remained strongly negative at the end of the measurement period.

Discussion
Our analysis demonstrates sustained long-term subsidence across a large 12 year dataset and supports previous findings that plantation management has contributed to subsidence both within the plantations themselves and in adjacent native forest (Hooijer et al., 2012;Evans et al., 2019).Multivariate analysis of the time series data suggests that water table depth exerts a dominant control on subsidence rates within

Table 4
Regression fits (Equations 5 and 6) to subsidence data for Acacia plantation, forest edge and forest interior sites, as shown in Figure 5. Regression fits to median Z score values for each group are also converted to annual subsidence rates for illustrative years.Total subsidence is calculated as the modelled change from 2007 to 2050 based on the regression equation shown, and does not include primary subsidence during initial forest to plantation conversion.both plantation and native forest areas.However, plantation sites were found to have slightly lower subsidence rates than forest sites (0.8 cm yr − 1 lower) if all other conditions were equal, and rates within plantations were slightly higher in areas closer to canals.These results suggest some additional influence of plantation management activities on subsidence, such as prior compaction of the peat during site preparation and harvesting operations, although these effects appear to be secondary.As noted earlier, we did not have detailed peat characterisation data for all monitoring sites so could not investigate this in detail, however it is likely that variations in peat type and properties account for some of the unexplained variance in the multivariate analysis.Time since initial plantation establishment was also found to exert a significant negative influence on subsidence rates in both plantation and native forest rates; this is discussed below.We found a clear influence of climate variations on subsidence rates.For the full dataset, higher rainfall was consistently associated with lower rates of subsidence.In the context of managed plantations, and to some extent also adjacent forest, rainfall and WTD may be partly decoupled due to plantation water management activities.Differences in the time period represented by each measurement (quarterly rainfall totals versus instantaneous WTD data) likely lead to further decoupling, although RAIN and WTD were nevertheless still correlated in the full dataset (r = 0.37), and it seems likely that most climatic impacts on subsidence are ultimately determined by their effect on WTD over weekly to monthly timescales.Based on the subset of sites with 12 years of monitoring data, it seems that large-scale climate indices provide a better representation of subsidence drivers than quarterly rainfall data, perhaps because they better capture periods of intense drought.In particular, the acute 2015-16 El Niño exerted a strong influence on subsidence rates, triggering an abrupt fall in peat surface elevation of around 4 to 10 cm (depending on land cover type) within a three month period.If we remove the effect of this short-term climate perturbation on the dataset by fitting regression models to the pre-2015 period, the estimated rate of subsidence in the plantation sites is reduced by 14 % (from − 5.17 cm yr − 1 to − 4.46 cm yr − 1 ).The compounding influence of the 2015-16 El Niño appears to have been smaller in the forest edge sites (7 % lower subsidence rates after removing the El Niño effect) but larger in the forest interior sites (32 % lower).
During the three years following the 2015-16 El Niño, it appears that the peat surface had largely rebounded to the level that would have occurred due to drainage alone in the plantation and forest edge sites; this is illustrated by the near-zero residuals at the end of the measurement period in Fig. 8a-b, and the convergence of the black (observed) and dashed red (linear model fitted to pre-2015 data) subsidence Z scores in Fig. 6a-b.On this basis, we infer that the somewhat lower estimates of long-term subsidence obtained by fitting empirical models to the pre-2015 data may provide a better estimate of the true rate of drainage-related subsidence than those obtained from a linear model fitted to the full dataset.
In the forest interior sites, the peat surface had not rebounded from the 2015-16 perturbation by the end of the record, as shown by the persistence of large negative residuals in Fig. 8c.This finding has a number of implications.Firstly, it suggests that the observed wholeperiod subsidence at these sites is at least partly associated with climate perturbation rather than with plantation impacts.This is consistent with the multivariate analysis, which did not show a clear influence of distance from the nearest canal on forest subsidence rates, and with a previous analysis of spatial variations in subsidence rates across the larger APRIL dataset (Evans et al., 2019) which did not show clear evidence of plantation impacts on water table depth or subsidence beyond around 300 m from the plantation edge.Secondly, our results suggest that native peat swamp forests are highly susceptible to severe and sustained subsidence during drought events.This may be explained by their typically lower bulk density compared to plantations (Anshari et al., 2010;Junedi et al., 2017), leading to rapid compaction and consolidation following loss of porewater from the upper peat layer.To some extent, peatlands naturally adjust to hydrological variations via oscillations in their surface elevation (sometimes termed 'bog breathing'), for example on a seasonal basis (e.g.Howie and Hebda, 2018).This self-regulatory function helps to reduce variations in peat surface wetness, conferring some resilience against climatic extremes.However, the magnitude of the elevation change in this case appears very large, and the extent to which observed ENSO/IOD related subsidence is reversible remains uncertain.Our results, to the end of 2018, suggest only slight recovery towards pre-2015 levels, and the El Niño event also appears to have a caused a sustained change in the relationship between peat surface elevation and water table (Fig. 7c).This could indicate that peat compaction due to the El Niño has led to an irreversible changeor at least very slow recoveryin the hydrological functioning of the peat at these sites.
At this stage, it is difficult to attribute observed subsidence in the forest interior sites to natural or human causes, either direct or indirect.It is possible, for example, that peat swamp forests are subject to periodic natural downward adjustments in surface elevation during dry periods, which are offset by the accumulation of low-density peat during wetter periods.This would lead to a long-term 'sawtooth' pattern of slow peat growth interspersed with short periods of more rapid subsidence.Even if this were the case, however, the amplification of ENSO/IOD impacts under a warming climate (Rifai et al., 2019) could make these periodic downward adjustments more frequent or severe.At worst, this could negate or reverse natural processes of peat formation, leading to degradation of peat swamp forests in the region even in the absence of direct human intervention.This interpretation also appears to be supported by recent eddy covariance CO 2 flux data collected over peat swamp forest in the Kampar Peninsular (Deshmukh et al., 2021) and in the Maludam National Park, Sarawak (Tang et al., 2020).These results, together with our subsidence data, could provide an early indication of the large-scale destabilisation of peat swamp forests due to intensifying climate change impacts, although more data would be needed to establish this.
With regard to direct human impacts, our previous observations of a limited (up to 300 m) lateral extent of water table drawdown into forests adjacent to plantations argue against a direct drainage-related contribution to subsidence in the forest interior (Evans et al., 2019).However we cannot rule out the possibility of more regional-scale impacts of landuse on the hydrology of peat swamp forests, for example by increasing topographic gradients at the edge of peat domes (Cobb et al., 2017) and subsequently reducing equilibrium peat dome volume (Cobb et al., 2020).The Kampar area included in our study (and that of Deshmukh et al.) may have been subject to hydrological modification by historic logging canals, while the study site of Tang et al. (2020) appears to have been affected by logging during the study period.Nevertheless, it seems clear that a large fraction of the observed subsidence at forest interior sites during our study period is not associated with steady subsidence due to plantation-induced water table drawdown.
With regard to the more sustained subsidence observed in plantation and forest edge locations, a key uncertainty for future plantation management, as well as for future CO 2 emissions (associated with the oxidative component of subsidence) and potential flood risks, is the trajectory of future subsidence.As noted earlier, long-term observations from drained temperate peatlands generally show exponentially decreasing subsidence rates over time, but no such records exist for tropical peatlands.Previous authors have argued that subsidence in these systems should proceed more or less linearly until most of the peat is lost (Hooijer et al., 2012;Couwenberg and Hooijer, 2013), whereas Hoyt et al. (2020) and Umarhadi et al. (2022) used satellite (InSAR) data and a space-for-time substitution to infer that subsidence rates were lower in areas that had been subject to longer-term drainage.Our dataset, which represents the longest continuous time series of subsidence observations from tropical peatlands, does indeed suggest possible slowing of subsidence over time, at least in plantation areas.Time since the start of monitoring was found to be a significant predictor of subsidence for all three land classes, consistent with sites that have been subject to drainage effects for longer having lower average subsidence rates (Table 1).For the subset of long time series plantation data, the R 2 of an exponential model fitted to pre-2015 data was marginally higher than the linear model (Table 2), again suggesting slowing rates of subsidence with time since drainage.This apparent slowing was however superseded by the rapid subsidence that occurred in response to the 2015 ENSO/IOD event, making any long-term projections uncertain.The selection of a linear versus a non-linear model made small difference to modelled subsidence during the short-term period, however increasing divergence would occur if the models were projected into the future.

Conclusions
Our analysis of one of the largest and longest-running subsidence monitoring datasets available for tropical peatlands reveals a profound influence of climate fluctuations on subsidence rates within both plantation and forest landscapes.In particular the severe 2015 dry season led to a sharp acceleration in subsidence rates, accounting for an estimated 14 % of overall subsidence observed over a 12 year period at plantation sites, rising to 32 % at forest interior sites.Whereas plantation and forest edge sites had, by the end of 2018, largely rebounded from the effects of the El Niño/IOD event, interior forest sites had not, raising the possibility that more frequent climate extremes could be contributing to slowly reversible or even irreversible changes in these ecosystems.Given the strong influence of the El Niño and IOD, determining the long-term trajectory of subsidence in plantation areas is difficult, however our analysis provides indications that subsidence rates is slowing with time since drainage, consistent with longer records from temperate peatlands.Given the low-lying nature of Southeast Asian peat landscapes, the magnitude of future subsidence will determine impacts on drainability and susceptibility to flooding.Our results suggest that the more extreme scenarios based on sustained linear subsidence may overestimate the long-term trajectory of peat subsidence.Nevertheless continued future subsidence could have major implications for future plantation water management, greenhouse gas emissions, and regional livelihoods.Resolving uncertainties in future subsidence projections for these economically important but vulnerable, carbon rich ecosystems requires continued and expanded subsidence monitoring, and an improved understanding of the interacting effects of plantation management and intensifying climate perturbations on peatland function.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 3 .
Fig. 3. Quarterly values of a) mean El Niño Southern Oscillation Index, b) Indian Ocean Dipole Index, c) quarterly rainfall totals averaged for all rain gauges on APRIL concessions, and d) median water table Z scores.The x-axis for the rainfall plot intersects the mean quarterly rainfall value for the measurement period.Median water table Z scores are based on all plantation and forest sites.

Fig. 4 .
Fig. 4. Quarterly median Z scores of water table depth for plantation (n = 38), forest edge (n = 13) and forest interior sites (n = 11) (based on the subset of sites with 12 years of monitoring data).

Fig. 5 .
Fig. 5. Median, 10th and 90th percentile peat surface elevation change relative to a December 2007 zero datum for Acacia, forest edge and forest interior subsidence monitoring sites.Gap in data for June 2014 at the Acacia plantation sites is due to missing data from over half of all monitoring sites at this time.Fig. 6.Median (solid black) and 10th/90th percentile (grey) standardised peat surface elevation data (Z-scores) for Acacia, forest edge and forest interior subsidence monitoring sites.Dashed lines show alternative regression models: Blacklinear fit to full dataset; Redlinear fit to 2007-2014 (pre 2015 El Niño) data; Greenexponential fit to 2007-2014 data.

Fig. 7 .
Fig. 7. Residuals of observed versus predicted quarterly subsidence Z scores obtained from a linear regression model fitted to measurements made prior to the 2015 El Niño, plotted against Z scores of water table depth, for each land cover category.Grey symbols show measurements made prior to the 2015 El Niño dry season, black symbols show measurements made after this time.Arrows denote an apparent step change decrease in peat elevation between September and December 2015.

Fig. 8 .
Fig. 8. Time series of quarterly residuals of observed versus predicted quarterly subsidence Z scores obtained from a linear regression model fitted to measurements made prior to the 2015-16 El Niño event.

Table 1
List of Variables for Multivariate Regression Analysis.

Table 2
Descriptive statistics of dynamic variables before and after outlier exclusion.