Short‐ and long‐term carbon emissions from oil palm plantations converted from logged tropical peat swamp forest

Need for regional economic development and global demand for agro‐industrial commodities have resulted in large‐scale conversion of forested landscapes to industrial agriculture across South East Asia. However, net emissions of CO2 from tropical peatland conversions may be significant and remain poorly quantified, resulting in controversy around the magnitude of carbon release following conversion. Here we present long‐term, whole ecosystem monitoring of carbon exchange from two oil palm plantations on converted tropical peat swamp forest. Our sites compare a newly converted oil palm plantation (OPnew) to a mature oil palm plantation (OPmature) and combine them in the context of existing emission factors. Mean annual net emission (NEE) of CO2 measured at OPnew during the conversion period (137.8 Mg CO2 ha−1 year−1) was an order of magnitude lower during the measurement period at OPmature (17.5 Mg CO2 ha−1 year−1). However, mean water table depth (WTD) was shallower (0.26 m) than a typical drainage target of 0.6 m suggesting our emissions may be a conservative estimate for mature plantations, mean WTD at OPnew was more typical at 0.54 m. Reductions in net emissions were primarily driven by increasing biomass accumulation into highly productive palms. Further analysis suggested annual peat carbon losses of 24.9 Mg CO2‐C ha−1 year−1 over the first 6 years, lower than previous estimates for this early period from subsidence studies, losses reduced to 12.8 Mg CO2‐C ha−1 year−1 in the later, mature phase. Despite reductions in NEE and carbon loss over time, the system remained a large net source of carbon to the atmosphere after 12 years with the remaining 8 years of a typical plantation's rotation unlikely to recoup losses. These results emphasize the need for effective protection of tropical peatlands globally and strengthening of legislative enforcement where moratoria on peatland conversion already exist.


| INTRODUC TI ON
The need for economic development across South East Asia, and global demand for agro-industrial commodities such as palm oil, rubber and pulp wood have driven the expansion of industrialscale agriculture and associated land-use change in recent decades.
Agricultural crop production now covers 122 million hectares in the region (Kenney-Lazar & Ishikawa, 2019), around a quarter of the total land area.
In Malaysia and Indonesia alone, oil palm plantations now cover an estimated 23 million hectares (Cheng et al., 2018;Gaveau et al., 2018;Miettinen et al., 2017). A significant proportion of this conversion has occurred recently on tropical peatlands; between 1990 and 2015, some 7.8 million hectares of these wetland peat swamp forests (PSFs) were converted through forest clearance and land drainage (Miettinen et al., 2016). The economic contribution of expanding oil palm production, particularly in rural areas (Qaim et al., 2020), has come at a, yet to be fully determined, cost to the local environment and the global carbon balance.
Measurements from soil surface chambers and subsidence focus on SOC losses which is often advantageous, but these measurements do not allow direct monitoring of the net CO 2 change in the atmospheric carbon pool as they cannot concurrently capture photosynthetic carbon uptake and above-ground sources of CO 2 . In this regard, the eddy covariance (EC) technique (Baldocchi, 2003) provides distinct advantages: EC measures the net ecosystem exchange (NEE) of carbon, capturing both emission and uptake, and spatially integrates over complex intra-site sources, such as drainage ditch peat extraction and autotrophic respiration from plant biomass. EC has previously been employed in peatland forests in the region; Hirano et al. (2012) used it to investigate the impact of large-scale anthropogenic disturbance on the carbon balance of tropical PSF in Indonesian Borneo, concluding that PSF were all now likely to be sources of atmospheric carbon (in the range of 7-18 Mg CO 2 ha −1 year −1 ). A conclusion supported by another, very recent, EC study in logged PSF that showed a mean net emission of CO 2 over 4 years at 15.4 Mg CO 2 ha −1 year −1 (Tang et al., 2020). However, the logistical and financial costs associated with EC have, to date, limited its deployment (Hill et al., 2017) and more studies are needed in tropical peatlands.
Only very recently have studies using EC started to report net carbon flux from oil palm plantations (Meijide et al., 2020), with only one monitoring oil palm cultivation on tropical peat (Kiew et al., 2020). The Kiew et al. (2020) (Mattsson et al., 2000;Schmidt, 2015).
In the absence of field studies of oil palm peatland conversions across the entire cultivation lifetime, emission factors determined for tropical forest conversion to agriculture on peatland have so far had to rely on very limited data. In deriving their Tier 1 emission factor of 40 Mg CO 2 ha −1 year −1 for conversions to oil palm on drained peatland, the IPCC list only eight direct studies, of which six were soil flux chamber studies and two were based on subsidence measurements. No ecosystem-level monitoring of carbon flux was available to be included in the assessment. The IPCC (Hiraishi et al., 2014) stated that emissions during the early years of plantation establishment are expected to be significantly higher than their emission factor but were not included due to this lack of available data.
The absence of directly measured net carbon flux from peatland conversions to industrial plantation led to controversy around the GWP impacts of conversion (Wijedasa et al., 2017). Reviewers for an Environmental Protection Agency (EPA) report into peatland emission factors for oil palm (EPA, 2014) were split over the importance of these early year emissions; debating the evidence of Hooijer et al. (2012) who had reported that very rapid subsidence recorded in the first 5 years of conversion was the result of large CO 2 emissions.
There was a suggestion that compaction may be contributing more to this large initial subsidence than the Hooijer et al. (2012) (Qaim et al., 2020), rely on two synthesis studies, Hooijer et al. (2010Hooijer et al. ( , 2012, as the key components in their assessment of the peat carbon impact of peatland conversion. The figure of CO 2 (Mg CO 2 ha −1 year −1 ) emissions being 91 times water table depth (WTD; m) from Hooijer et al. (2010) is the default value in their current GHG calculator (https://rspo.org/certi ficat ion/palmg hg/palm-ghg-calcu lator) which underpins estimates of peat decomposition in their certification scheme (https://rspo. org/certi fication). This coefficient of 91 was derived from a linear fit (with a fixed intercept of zero) to just eight data points collated from five studies, yet plays an important role in modelling the global carbon budget (Friedlingstein et al., 2019;Houghton & Nassikas, 2017;Le Quéré et al., 2018).
In this study, we begin to address this important knowledge gap by contributing data collected by EC at two adjacent oil palm plantations established on tropical peatlands in South East Asia.
One site captures a period following the initial conversion from PSF and another captures the mature phase. We present annual net ecosystem CO 2 fluxes from individual measurement years at both sites and partition them into photosynthetic uptake and whole ecosystem respiration. We then combine both data sets into a single chronosequence over a 151-month period and use a mass balance approach (incorporating estimates of biomass accumulation and forest residue decomposition) to calculate changes in soil carbon stocks. We present emission factors both for individual years and across relevant time periods (e.g. years 1-6, as highlighted by the IPCC). Finally, we investigate the relationship between soil water drainage and carbon loss in the context of previous emission coefficients and consider the potential impact of changes in plantation drainage targets.

| Site location and description
The two study sites were individual blocks of commercially managed oil palm plantation situated within the Sabaju (OPnew:

| Instrumentation
Eddy covariance was carried out at both sites using identical instrumentation with the only significant difference being that profile measurements, for canopy storage of CO 2 and energy, were from three heights on a 20-m tower for the taller palms at OPmature compared to two on a 6-m tower at OPnew. LI-COR closed path systems (LI-7200/7550; LI-COR Environmental coupled to R3-50 Sonic Anemometer; Gill Instruments Ltd.) were used at both sites, with sensors sited at the top of each tower. For OPmature this resulted in a measurement height (above-ground) of 20.19 m, approximately 12 m above an 8 m canopy; for OPnew, sensors were at 6.06 m above a canopy that reached 2.6 m by the end of the study period. Prior to canopy development at OPnew topography was dominated by forest destruction residues compacted into rows of approximately 2 m in height which gave a typical measurement height above canopy of around 4 m. Canopy profile CO 2 and energy storage was measured using CO 2 diffusion sensors coupled with air and relative humidity sensors (GMP343 and HMP155A; Vaisala Corporation). For OPmature these were placed at 1, 6 and 18 m above-ground, for OPnew this was at 1 and 6 m. Energy balance was monitored at two locations for each site using heat flux plates (HFP01SC; Hukseflux Thermal Sensors) at 0.08 m soil depth coupled to soil moisture/temperature sensors (Steven's Hydraprobe; Stevens Water Monitoring Inc.) at 0.04 m. WTD was monitored within 0.05-m-diameter porous plastic pipe inserted to a depth of 2.5 m (PX709GW submersible pressure transducer; Omega Engineering Inc.). Precipitation was measured at the top of each tower using a tipping bucket gauge (TR-525M; Texas Electronics). EC data were collected at 10 Hz and written to an industrial-grade USB drive within the LI-7550, meteorological data at 1-minute intervals stored to Xlite 9210 dataloggers (Sutron Corporation).
Statistical outliers (spikes) in the 10 Hz data were detected following Vickers and Mahrt (1997); vertical wind speed measurements were only accepted at <5 SDs (σ) from the 30-min mean, other variables at 3.5σ; 30-min periods containing spikes at greater than 1% were flagged as poor quality. Time lags, discrepancies between precise sampling times at the anemometer and gas analyser, are compensated for using site-specific covariance maximization derived from data collected at the site. Detrending of turbulence fluctuations over each 30 min was through block averaging. Co-ordinate rotation, to accommodate imperfect alignment to the horizontal wind vector, was carried out through the planar fit method of Wilczak et al. (2001) using sitederived parameters. Co-spectral analysis and correction of low-and high-pass filtering effects were carried out following Moncrieff et al. (1997Moncrieff et al. ( , 2004. Spatial estimation of the areal source of sensor data capture (footprint assessment) followed Kljun et al. (2004) or Kormann and Meixner (2001), where turbulent friction was <0.2 m s −1 . CO 2 storage below the sensor height, which is not captured in turbulent eddy transfer through the EC sensor pathways, was accounted for through profile monitoring of CO 2 concentrations between the ground and sensor heights as outlined in Baldocchi et al. (2001): time-stamped changes in absolute CO 2 concentration are captured by profile sensors and converted to volumetric ratios using the ideal gas law. These are then added to the corresponding flux measurements captured by EC at the half-hour time step.

| Quality control flagging
Initial quality control flagging of each 30-minute flux average (statistical testing of the 10 Hz data) followed the Carbo-Europe standard 0-1-2 system of Foken et al. (2004). Zero being the highest quality, values flagged at 2 were automatically discounted from further processing. Data spikes in the half-hourly processed CO 2 data were identified and removed following Papale et al. (2006) using the sug-

| Study site area (measurement fetch and footprint)
The available study area which satisfied EC assumptions of homogeneity and representation of the area of interest (fetch) covered 41.7 ha at OPnew and 907 ha at OPmature. A combination of Google Earth (GE v7.3.2.5776, Imagery date 24/03/2016) and ARCGIS (ArcMap 10.5.1; ESRI) was used in conjunction with the output from the footprint model to filter out any measurement periods where data collection extended beyond the ideal fetch. Taking the sensor tower location as a datum point, distances to the edge of the fetch boundary were measured at 10° increments. Half-hourly output from the footprint model (percentage data contribution to total, distance to peak contribution and wind vector) was compared to these boundaries within 10° bins (total of 36) and considered acceptable where 70% of the information collected in each half-hour was sourced within the fetch boundary.

| Energy balance
Energy balance closure (EBC) was investigated using an ordinary least square (OLS) regression at the half-hour time step between turbulent heat flux (LE plus H) and available energy (net radiation plus soil heat flux). EBC was considered as the slope of the resulting OLS fit. Energy storage in relevant pools (canopy air space and soil volume) was calculated using specific heat capacity and moisture fluctuations. Energy lost to photosynthetic utilization was calculated following Masseroni et al. (2014). The ratio between turbulent heat flux and available energy over the entire study period is presented as the energy balance ratio following Wilson et al. (2002).

| Gapfilling and flux partitioning
Gapfilling of data rejected through quality control and partitioning of NEE into photosynthetic uptake (gross primary productivity [GPP]) and ecosystem respiration (R eco ) was carried out using the ReddyProc package (Wutzler et al., 2018) within R. For gapfilling, this package utilizes the mean diurnal separation (MDS) approach of Falge et al. (2001) with flux partitioning carried out using the light response curve method of Lasslop et al. (2010) to estimate daytime GPP, the sum of NEE and GPP being R eco . Night-time fluxes (below a global radiation (R g ) threshold of 20 W m −2 ) are assumed solely R eco (Reichstein et al., 2005). Underestimation of fluxes during periods of insufficient turbulence was avoided by removing data recorded below site-derived friction velocity thresholds (u* filtering) during the gap-filling process (Reichstein et al., 2005). The first 4 months of data immediately following the conversion of PSF at OPnew were not collected due to the sensor installations not yet being in place. An estimation of these values has been made through modelling backwards at a monthly time step from the first data available. For GPP, an assumption is made that this would be zero immediately following conversion (dead forest residues and bare soil), therefore linear interpolation was carried out from a start point of zero at the beginning of May 2016 to the beginning of September 2016 (the first complete month's data). For R eco , a linear trend line was fitted to the existing R eco monthly data set and extended back over these first 4 months.
Trends in fluxes (GPP, R eco and NEE) over the measurement period at each site are indicated by the slope of a linear model fitted against time, with significance accepted at p < 0.05. Estimation of annual means within specific periods is calculated by multiplying the mean monthly values for that period by 12. A complete chronosequence of NEE, from months 1-151 was established using exponential interpolation (Stineman, 1980) between the OPnew and OPmature data sets at a monthly time step. This interpolation was also applied to months 141-146 (April 2019-August 2019) at OPmature, where data were excluded due to sensor malfunction resulting in a gap too large for the MDS gap-filling routine.

| Live standing biomass
Net primary productivity (NPP) is the sum of photosynthetic carbon sequestered into biomass pools on site during any given period. For live palm biomass carbon stocks, data were interpolated at a monthly time step between biomass, and associated carbon concentration, for age classes 3, 8 and 12 years, measured from destructive sampling from appropriately aged planting blocks at the Sabaju (age classes 3 and 8 years) and Sebungan plantations (age class 12 years) in 2019 and presented in Lewis et al. (2020).
Individual palm component carbon stocks are summed to total palm biomass carbon using Equation (1). Root biomass was not directly sampled in Lewis et al. (2020) so has been assumed at 16% of total standing biomass following Khalid et al. (1999). This where all components' dry mass multiplied by fraction of carbon content (see Table S

| Fresh fruit bunch harvest offtake
Harvest offtake, as FFB, was provided by the site managers for the OPmature planting block specifically within the Sebungan plantation at a monthly time step from the date of first harvest in month 32 to month 144. For months 0-31, linear interpolation was used to complete the monthly timeseries from zero to first harvest. For sequestration calculations, all FFB is considered to remain within the system (see Section 4), therefore total NPP of FFB for any given period is the cumulative sum of all harvest offtake during that period.

| Pruned frond biomass
At each harvest, a number of fronds are cut to facilitate access to FFB and left in piles to decompose on site. While uptake of carbon into these fronds during growth, and return to atmosphere through decomposition, will be captured by EC in NEE, the carbon stored within the ecosystem in the total frond pile biomass at any given time needs to be accounted for in NPP. An estimate for this was derived from monthly interpolation between the numbers of frond bases (remnants of removed fronds) present on the palms at the 3-, 8-and 12-year time points. While multiplying these by a mean frond mass (for each age class of frond, similarly interpolated) gives an estimate of frond mass pruned in each month, account needs to be taken of the decomposition of each monthly addition to the pruned frond pile over the remaining study period. This was carried out by applying an exponential decay function, Equation (3) (Moradi et al., 2014;Olson, 1963) to each monthly pruned frond mass and continuing to the end of the chronosequence, then summing across the remaining biomass from all previous prunings to a total pruned frond biomass pool per hectare for each month. The decomposition rate constant (fractional mass loss per month [k]) for frond biomass was set at 0.15, calculated from empirical measurements by Moradi et al. (2014).
where t denotes time (month) and k denotes decomposition rate constant (fractional mass loss per month).

| Calculation of forest debris decomposition
The contribution to ecosystem respiration (R eco ) from the decomposition of the previous forest biomass (R fr ) that was cut and compacted on site prior to establishment of oil palm needs to be considered in the overall carbon budget as these emissions will be a significant contribution to the net flux captured by EC. Starting biomass and decomposition rate were not measured directly on site, so literature estimates have been used. As with the frond pile biomass, the decay rate for forest coarse woody debris (CWD) decomposition is calculated using Equation (3) at a monthly time step, but in this case just considers a single biomass addition at the beginning of the conversion. Kho and Jepsen (2015) estimated 58.7 ± 10.7 Mg C ha −1 for logged PSF (as at OPnew); using their dry stem biomass carbon content of 47.4%, we derive forest biomass at clearance of 123.8 ± 22.6 Mg DM ha −1 . Only one paper was found that monitored CWD decay under tropical peatland conditions in the same region (Law et al., 2019), who were working in Sabah, Malaysian Borneo under very similar climatic conditions to this present study. That study reported that after 12 months there had been a mean loss of (1) Carbon biomass = Carbon roots + Carbon trunks + Carbon frondbases + Carbon fronds +Carbon FFB + Carbon spears + Carbon cabbage + Carbon residual , 25.6% of the starting biomass. From this we calculated a rate constant (k) at 0.02464 for decomposition of forest CWD at a monthly time step which would result in the correct biomass loss by month 12.

| Changes in SOC
Changes in SOC (ΔSOC) during specific periods were calculated as the difference between GPP and R eco , taking account of sequestration of GPP into biomass stocks (NPP) and contribution of R fr to R eco . Following production of a complete monthly timeseries of NPP and R fr , as described above, periods corresponding to direct measurements of NEE (partitioned into GPP and R eco ) were used as parameters in Equation (4) to estimate peat carbon loss for the study period months at OPnew and OPmature. For annual carbon emission factors over the entire 151-month chronosequence of NEE, and for specified periods, Equation (4) is modified using Equation (5) (see Equation S.2.1) and uses parameters derived by taking mean monthly values within selected time periods and multiplying by 12.
Emission factors were calculated for years 0-6 (establishment), 6-12 (mature) and across the entire period. Additional exports of carbon (ε), such as drainage losses of dissolved and particulate organic carbon or methane emissions (CH 4 ), are not captured in our EC results and are therefore not accounted for in these calculations; however, their potential magnitude is considered in the discussion section below.
where ΔSOC is the change in soil organic carbon (Mg CO 2 -C ha −1 ), GPP is the gross primary productivity (photosynthetic uptake of carbon, Mg CO 2 -C ha −1 ), R eco is the Ecosystem respiration (Mg CO 2 -C ha −1 ), NPP is the net primary productivity (carbon sequestered into biomass, Mg CO 2 -C ha −1 ), R fr is the respiration contribution from decomposition of forest residue (Mg CO 2 -C ha −1 ), ε denotes the unaccounted factors (e.g. export as dissolved organic carbon, carbon content of emitted CH 4 , etc., Mg CO 2 -C ha −1 ) and t denotes the time (year) where NEE is the net ecosystem exchange of carbon (Mg CO 2 -C ha −1 )

| Relationship between peat carbon loss and WTD
The relationship between our estimate of SOC emission (as CO 2 ) and WTD was considered in two separate analyses.  2006)) and subsequent inclusion of annual ΔSOC from OPnew and OPmature. Secondly, taking advantage of our high-frequency respiration data (at the half-hour time step) at OPnew, we fit a non-linear second-order polynomial curve to nighttime NEE, assumed to be entirely respiration, and WTD following detrending of the data series and binning of R eco into 0.01 m increments of WTD. Only original measured data (not gapfilled) were used and selected at the highest quality (qc flagged at zero ;Foken et al., 2004). The output of the model fit was then used to predict respiration rates for 0.1 m WTD increments between 0 and 0.8 m below the soil surface (the measured WTD range at OPnew within the study period).

| RE SULTS
See Supporting Information for full details of EC data capture and retention following quality control. Values ± given throughout these results indicate the standard error of the mean (SEM)

| Climate data
As expected in equatorial, tropical climate rainfall was high and temperatures were relatively stable over time with only a minimal seasonal component (Figure 1). Rainfall averaged 2856 ± 96 mm year −1 across the two sites, and had a mean temperature of 26.9 ± 0.03°C.
The exposed soils at OPnew were on average 2°C warmer than at

| Soil organic carbon loss (as CO 2 -C) for individual study years at both sites
The ΔSOC was 2.5 to three times higher in OPnew compared to

| Chronosequence of cumulative NEE
Total cumulative NEE across the entire chronosequence suggested a net emission of CO 2 from the site at 823.3 ± 0.9 Mg CO 2 ha −1 (224.9 ± 0.2 Mg CO 2 -C ha −1 ) after 151 months. As can be seen in  Abbreviations: GPP, gross primary productivity; NPP, net primary productivity; R eco , total ecosystem respiration, NEE, net ecosystem exchange of carbon; R fr , respiration of carbon from decomposition of forest residue; ΔSOC, changes in soil organic carbon during period.  G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G  Cumulative NEE (Mg CO 2 ha −1 )

| Annual carbon emission factors (as CO 2 -C)
Mean annual NEE across the 151-month chronosequence showed a net annual emission of 17.9 ± 1.3 Mg CO 2 -C ha −1 year −1 (65.6 ± 4. 8 Mg CO 2 ha −1 year −1 ). Using Equation (4) (modified by Equation 5) to account for carbon sequestered into on site biomass at a mean annual NPP of 4.9 ± 0.2 Mg CO 2 -C ha −1 year −1 and a mean annual carbon emission from the decomposition of forest residue (R fr ) at 4.6 ± 0.1 Mg CO 2 -C ha −1 suggested a mean ΔSOC over the entire period at −18.3 ± 1.3 Mg CO 2 -C ha −1 year −1 (equivalent to a soil surface emission of 67 ± 4.8 Mg CO 2 ha −1 year −1 ).

| Linear fit of ΔSOC (as CO 2 flux) to WTD
Adding our sites' annual soil carbon emissions (ΔSOC in Table 1, as CO 2 ) and mean WTD (see Figure 4)

| Non-linear fit of ecosystem respiration to WTD
The second-order polynomial fit of night-time R eco to WTD (see  Table 3 shows predicted annual R eco for each 0.1 m WTD increment and their magnitude relative to a typical plantation drainage target of 0.6 m. Figure 5 shows a graphical representation of these estimated percentage changes in annual R eco when comparing different mean annual WTDs. For example, raising WTD to 0.2 m would see a 31% reduction in respiration compared to the typical 0.6 m target. In  contrast, differences between depth increments deeper than 0.6 m were insignificant.

| CO 2 flux
We have presented, for the first time, a comparative study of measured NEE of carbon as CO 2 between the initial years of peatland conversion to oil palm and the later mature phase. Our results show the dramatic difference between the very high early soil carbon emission remained high at around 47 Mg CO 2 ha −1 year −1 but was being masked by NPP. Sequestration of carbon into the biomass pool (NPP) was equivalent to around 54% of the peat carbon loss for that period. A small proportion of that carbon will be contained within the fruit harvest offtake, removed from the site and returned to the atmosphere during oil production and consumption, the remainder will be held within the palm and frond litter biomass pool until re-cultivation (typically at around year 20) where it will begin to return to the atmosphere during the decomposition of the palm biomass following clearance.
We also considered the relative contributions to NEE from uptake (GPP) and emission (R eco ). The common approach for partitioning NEE (Reichstein et al., 2005) uses parameterization of the relationship between night-time NEE (assumed to be entirely R eco ) and air temperature to estimate the contribution of R eco to daytime NEE (the residual being GPP). However, this approach relies on a strong relationship between respiration and temperature. This can be problematic in tropical climates where temperature ranges (both diurnally and seasonally) tend to be very much narrower than in temperate zones, and likely compounded by the strong relationship in drained peatlands between R eco and WTD. As an alternative, we adopted the approach of Lasslop et al. (2010) using a light response curve fitted to daytime NEE to estimate GPP. Kiew et al. (2020) concluded that low GPP in their poorly established plantation was responsible for their large on site net emissions. In line with their conclusion, as seen in Figure 2, while R eco was slightly lower at our mature site compared to OPnew, it was the much higher GPP into the mature palms that was primarily responsible for driving this reduction in NEE. There was a small, but statistically significant, reduction in both GPP and R eco over time during the OPmature monitoring period, which appears to be most apparent in the period between

| Relationship between WTD and soil carbon loss (as CO 2 )
Adding our estimate of soil carbon loss (as CO 2 ) from OPmature to the Hooijer et al. (2010) linear regression model, we found that our mature site fitted remarkably well within their original data set, only raising the coefficient to 93.7 from their CO 2 = 91*WTD. However, including our early conversion period at OPnew increased this sensitivity estimate by 26%. This reinforces the importance of incorporating the early years of conversion into assessments of the carbon impacts of peatland conversion and emphasizes the need for emission factors covering entire cropping periods.
A single coefficient such as this may, though, be too simplistic. As discussed in Hooijer et al. (2006), a coefficient for the relationship between CO 2 flux and WTD is unlikely to be consistent throughout the soil profile, and it was a lack of available data which limited their original analysis to a simple linear fit. In their discussion (see note under figure 12 in Hooijer et al. (2006)), the authors suggest that CO 2 emissions might be reduced at WTD of 0.2-0.3 m and at zero when WTD = 0, that is, waterlogged conditions would promote the formation of peat and net CO 2 emissions would <= 0. Their suggestion that a linear coefficient of 91 would hold true for WTD from 0.25 and 1.1 m was discussed in Couwenberg et al., 2010) who concluded that there was not enough evidence to clearly state whether rates of subsidence (as a measure of peat decomposition) became static beyond a WTD of 0.5 m.
Taking advantage of our high-frequency, long-term data set at OPnew, where heterotrophic respiration should heavily outweigh the limited autotrophic respiration from immature, widely spaced palms, we investigated this relationship using a nonlinear, polynomial fit. Here we found not only a highly significant relationship (in contrast to Cooper et al., 2020, whose study had very little regression data available) but one that implied the opposite to what was suggested in Hooijer et al. (2006). We found a much greater sensitivity to drainage in the upper half of the soil profile compared to the lower half over our 0-0.8 m WTD range (see Table 3). Each additional 0.1 m drainage within the upper 0.4 m produced an additional 24.5 Mg CO 2 ha −1 year −1 (compared to 9.1 Mg CO 2 ha −1 year −1 expected from Hooijer et al. (2010)), while between 0.5 and 0.8 m this decreased to 6.8 Mg CO 2 ha −1 year −1 . This trend may be reasonably intuitive, Leifeld et al. (2012) showed that peat decomposition rates were dependent on organic matter quality and that decomposability was higher in the newer organic matter nearer the peat surface. Given the importance of temperature in driving soil carbon decomposition (Lloyd & Taylor, 1994), it would also follow that However, more work is needed to investigate the impact that reducing WTD would have on fruit yield. Hooijer et al. (2012) estimated (from subsidence studies) that the first 5 years of conversion from PSF to plantation (acacia as this is the conversion they used to estimate 0-5 year fluxes) would see a mean loss of peat carbon at 48.6 Mg CO 2 -C ha −1 year −1 (calculated from CO 2eq .

| Comparison to subsidence studies
figures in their Table 2 It should also be noted that our estimate of CO 2 flux from the decomposition of the forest biomass was contributing around a quarter of the total ecosystem CO 2 emission to the atmosphere (R eco ). This contribution would be reducing the estimate of early years' peat decomposition considerably in our mass balance equation, 75% of the forest residue was decomposed within 4 years in our decomposition model. This estimate was based on a literature figure for the starting biomass (Kho & Jepsen, 2015) and an assumption that 25% had decomposed by the end of the first year from a single decomposition study (Law et al., 2019 An aspect not captured in our direct monitoring of the ecosystem/atmosphere exchange of CO 2 is peat loss due to the export of carbon in groundwater as dissolved and particulate organic carbon. This is represented within ε in Equation (4)

| Comparison to existing emission factors (as CO 2 )
Our long-term emission factor (calculated across all 151 months) for peat carbon loss (as CO 2 ) at 67 Mg CO 2 ha −1 year −1 is closer to the IPCC emission factor for peatland conversion to acacia plantation, 73 Mg CO 2 ha −1 year −1 than oil palm which is lower at 40 Mg CO 2 ha −1 year −1 (Hiraishi et al., 2014). The effect of plantation species on peatland CO 2 emission was not found to be a significant factor in the study of Carlson et al. (2015), who discussed the likely importance of time since drainage, though their data set was limited to a narrow age range. Miettinen et al. (2017) preferred to use the mean of the two IPCC factors (55 Mg CO 2 ha −1 year −1 ), in their calculation of carbon loss across the region due to peatland conversion which agrees well with the Cooper et al. (2020) mean figure of 53.1 Mg CO 2 ha −1 year −1 . All these estimates remain lower than the EPA-accepted emission factor of 95 Mg CO 2 ha −1 year −1 (EPA, 2014). The IPCC explicitly exclude the first 6 years of conversion in their emission factor due to lack of data but acknowledge that this period would see much higher carbon losses. This is clearly demonstrated by our estimate of soil carbon emission for years 1-6 at 91.6 Mg CO 2 ha −1 year −1 . The overall net emission of CO 2 to the atmosphere (NEE) for this period, incorporating forest biomass decomposition and photosynthetic uptake, was higher at 110.8 Mg CO 2 ha −1 year −1 .
This difference between direct soil carbon emission (ΔSOC, as CO 2 ) and the net ecosystem scale addition of CO 2 to the atmosphere (NEE) is an important distinction, particularly in the early years of conversion, and should be considered when assessing the impacts of land-use change. In Table 2, we present, based on our results, estimates of emission factors for both NEE and ΔSOC at an annual time step across a 12-year period. Taking the cumulative sum (for either component) for any given period and dividing by that number of years will provide an estimate of mean annual CO 2 emission for that period.

| CON CLUS IONS
Despite our results reporting lower peat carbon loss in the early years following conversion than subsidence studies might have suggested, there is no doubt that these emissions remain extremely significant and PSF conversion to agriculture results in very large net emissions of CO 2 . We have shown that the impact of these fluxes on the atmospheric carbon pool can be larger than emission factors for soil carbon loss alone might suggest, and that is highly unlikely that 'carbon debts' incurred early in a plantation lifecycle could be recouped over the remaining years. Evidence has shown that moratoria on peatland conversion within protected areas can be effective (Chen et al., 2019) but huge challenges remain, and newly identified areas of extensive peatland being reported from tropical zones across the globe (Lähteenoja & Page, 2011;Xu et al., 2018) reveal regions that may be particularly vulnerable to land-use change. Despite policy development in South East Asia, limitations in regulatory frameworks and enforcement capabilities combined with political and socio-economic factors still challenge peatland protection (Padfield et al., 2016;Wijedasa et al., 2018).
Our results should make it clear that conservation of these globally important carbon stocks is vital to any efforts to minimize the impacts of future climate change and reduce the contribution of land-use change to it.

ACK N OWLED G EM ENTS
The authors thank the Director-General of the Malaysian Palm Oil Board (MPOB) for permission to publish these results. This study was carried out as part of a wider tropical peat research collaboration between MPOB, University of Exeter, University of Aberdeen and Newcastle University, and they thank the Sarawak Oil Palm Berhard (SOPB) for their help and support during the project.
Specifically, at SOPB they thank Paul Wong Hee Kwong (group CEO), Chua Kian Hong (group plantation manager), Phang Seng Nam (regional plantation controller) and the Sabaju and Sebungan plantation managers for being kind enough to allow the research platform to be established within their plantations and for providing logistical support when needed. At MPOB, they particularly thank the dedicated field technicians, without whose efforts and commitment the research would not have been possible, spe-

CO N FLI C T O F I NTE R E S T
The research was carried out as part of a project funded by the Malaysian Palm Oil Board (MPOB). L. K. K and E. R. are both employees of MPOB. The research was carried out with the support of Sarawak Oil Palm Berhard (SOPB) on whose land the research project was based.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.