The data were collected during August and September 2014 from the Sebungan oil palm plantation in Sarawak, Malaysia (latitude 003^o09’ N, longitude 113^o21’ E). Sebungan Estate has been established on 4.0 m deep peat soils broadly classified as histosols (FAO, 2006). The plantation has a tropical climate; the mean annual temperature was 26°C and the mean annual precipitation was approximately 3,000 mm (Cook et al., 2018; McCalmont et al., 2021). The northeast monsoon from October to January has the most rainfall, with a slightly drier southwest monsoon between May and August (Cook et al., 2018).

Prior to planting, the land use was a mixed species swamp forest, which had been heavily logged. The land was converted to a plantation in 2006 and the palms were on their first crop rotation. The palms were 7 years old when measurements began. The plantation was laid out systematically with ~35 ha blocks and drainage ditches every 28 m leading to a larger ditch running down the center of the block. Palms were planted every 8 m in rows that were 8 m apart, leading to a planting density of 160 palms per ha. Within the palm blocks, four different surface management microforms were present and two different drain types (Figure 1):

Six palms were selected within a 1 ha plot for this experiment. Three palms were located in rows next to the field drains and three palms were located in rows next to the rows of cover plants. At each palm, three sampling transects were set up, with each transect going across a different surface management microform (Figure 1). These transects within each management microform enabled us to determine if the effect of management microform interacted with root density to influence R_tot, R_h and R_a.

R_tot measurements were collected on the 30th and 31st August 2014. Samples were collected at 0, 0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 3.5, 4 m from the base of each palm. A further sample was taken at a point equidistant from the three nearest palms (i.e., 4.5 m from the palm bases). Flux measurements were performed with a static chamber approach, using 10 cm diameter flux chambers (Livingston and Hutchinson, 1995). Chamber bases were installed in the soil to a depth of 5 cm 4 weeks prior to the commencement of sampling in order to avoid disturbance effects associated with base installation. R_tot measurements were made in triplicate using a PP Systems EGM-4 and SRC-1 chamber (Hansatech Instruments Ltd., Norfolk, UK). For each replicate measurement, CO₂ concentrations were measured over a 2 minute enclosure period, with concentrations recorded at 3 second intervals, or until an increase of 50 ppm CO₂ had been observed.

Ambient air temperature was measured at the same time as the R_tot measurement using a thermometer (LCD Digital Thermometer, ATP Instrumentation, Leicestershire, UK; precision ±1°C). Soil temperature and soil moisture measurements were taken following the completion of R_tot measurement, adjacent to the collar as in Marthews et al. (2012) and Manning et al. (2019). Soil moisture was measured using an ML3 probe and HH2 moisture meter (Delta-T, Cambridge, UK; precision 1%). Following the chamber measurements for each palm, water table depth (WTD) was determined by digging a hole in the peat in the harvest path to the water table, 2 m away from the palm.

After the flux measurements, a 10 cm diameter, 30 cm deep, soil core was collected from each collar. Root dry mass was determined in each soil core using the method developed by Metcalfe et al. (2008), sampling up to 50 min for each soil core. The roots were washed and dried at 70°C to constant weight (Metcalfe et al., 2008). Root density was calculated as mass of dry roots over total volume of core extracted.

Following the extraction of roots, soil samples were dried and sieved to 2 mm. Chemical analyses were performed at the Malaysian Palm Oil Board headquarters in Kuala Lumpur, where soil C and N were measured in a CNS analyzer (Elementar Vario MACRO Cube, Germany), and soil pH was determined (Thermo Orion pH/ORP/cond model 555A, Thermo Fisher Scientific, Chelmsford, MA, USA).

Bulk density was estimated by finding the dry mass of soil with known volumes, sampled from the top 10 cm in the harvest path, frond piles and cover plants with bulk density rings. Soil C and N content were estimated by multiplying the bulk density by the soil C and N for the top 10 cm by 1 m of soil.

R_tot fluxes were calculated using R version 2.15.1 GUI 1.52¹. Here linear regressions were fitted and flux estimates determined (Equation 1) using the method in Marthews et al. (2012):

where ΔC is the change in CO₂ over the measurement period (ppm), P is pressure (mb), V is volume (m³), M is the relative molecular mass of CO₂, Y is the conversion to upscale the flux to annual emissions, Δt is the duration of the measurement period (s), T is temperature (K), A is surface area (m²) and R is the Universal Gas Constant 8.31432 J mol⁻¹ K⁻¹.

R_tot was partitioned into R_h and R_a using Equation 2:

Two different methodological approaches were used to estimate the partitioning of R_tot into R_h and R_a:

For the analysis, the aggregated mean R_tot per distance class per palm and aggregated mean root density per distance class per palm were used, reducing the maximum n = 210 to n = 72. This was done because it had been estimated that the proportional representation of each microform within each distance class ring was equal by creating a spatial model of the oil palm plantation and determining the size of each land cover type.

Estimating R_h and R_a at each individual measurement points was modeled using the linear regression in section 2.8. Here, measured root density was substituted into the linear regression equation in order to model R_a. Equation 2 was used to model R_h for each collar.

R_tot, R_h and R_a fluxes were upscaled to plantation-scale estimates using two different approaches: straight mean averaging and area-weighted upscaling.

The standard deviation was used to determine the confidence intervals (CIs) for the straight mean averaging estimate of R_tot and the straight mean averaging estimate of R_h and R_a calculated using the distance from palm method. It was possible to scale these CI for the distance from palm area-weighted upscaling estimate. Confidence intervals for R_h calculated using the linear regression method were found by estimating confidence intervals for the linear regression line and the estimates at the intercept were used for the R_h estimate. CIs for R_a using the linear regression method were found using the combination of errors technique.

These methods of determining CIs could not be applied to area-weighted upscaling estimates of R_tot or R_h and R_a estimated using the linear regression method. In these instances, CIs were determined by resampling the data in R 10,000 times, repeating the upscaling analysis and determining the 5th and 95th percentile mean estimate. CIs for the distance from palm area-weighted upscaling R_a were found using the combination of errors technique.

All statistical analyses were performed using R version 2.15.1 (see footnote 1). The estimates of R_h and R_a used for statistical analysis were taken from the method described in section 2.9 in order to obtain individual points.

Linear models were used to determine how R_tot, R_a and root density varied with distance from palm and with sampling transect. Tukey HSD tests were used for post hoc analyses, using the ANOVA function in R. Kruskal-Wallis tests were used for R_h. Nemenyi tests using Tukey HSD were used for Kruskal-Wallis post hoc analyses (R package: PMCMR v4.3; Pohlert, 2018). Following each statistical model, the fixed effect (and where necessary also random effect) residuals were considered for normality using Shapiro–Wilk normality tests, heteroscedasticity and equal variance.

Linear models were used to determine whether environmental variables were significantly different between distance from palm or the different surface management microforms. A log transformation was used for soil moisture and Box Cox transformations were used for soil carbon (lambda = 3.9), soil nitrogen (lambda = −2.05), soil C:N (lambda = 1.55) and soil pH (lambda = 3) to achieve normality with the model residuals.

The effects of environmental variation on R_tot were considered using a linear mixed effect model using the nlme package in R (Pinheiro et al., 2017). The fixed effects in the model included soil moisture, soil temperature, air temperature, root density and an interacting fixed effect factor between soil pH and soil surface microform transect. Distance from palm was included as a random effect.

Variation in the plantation-scale estimates of R_tot, R_h and R_a were assessed based on:

In order to assess variation in R_h due to partitioning methodology, the variation in R_h estimates produced in this study was considered. To assess sample size, subsets of the data were resampled from the dataset at sample sizes 5, 10, 20, 35 and the entire dataset. R_tot, R_h and R_a were then estimated using straight mean averaging. Furthermore, a power analysis was applied to the dataset to determine the number of samples needed to accurately estimate R_tot, using the method in Metcalfe et al. (2008). Here Equation 3 was applied to the dataset:

where t_α is the statistical significance wanted for the power analysis (here 0.05), CV is the sample coefficient of variation, and D is the specified confidence interval (here 10; Hammond and McCullagh, 1978).

With the intention of assessing variation in R_tot, R_h and R_a due to sampling design, random sampling with straight mean averaging and spatial sampling with area-weighted upscaling were compared. To show the variation in R_tot, R_h and R_a whether the rhizosphere was included or not, plantation-scale estimates of R_tot, R_h and R_a were estimated with and without rhizosphere data. Similarly, to consider whether including or excluding the different surface microform data, plantation-scale estimates of R_tot, R_h and R_a were considered for the different sampling scenarios with the harvest path data only or with the full dataset (Supplementary Tables S1, S2).

Best estimates of plantation R_tot, R_h and R_a are presented in Table 3. Straight mean averaging results were higher than estimates based on area-weighted upscaling. Distance from palm gave higher R_h estimates and lower R_a estimates than the linear regression method. Estimates of R_a gave greater changes with upscaling method than R_h estimates.

Methodological decisions influenced the estimates of R_tot, R_h and R_a. The sample size of data points impacted the final R_tot, R_h and R_a estimated from oil palm plantations on peat soil. A power analysis showed that 35 samples are needed to accurately capture the within-plantation spatial variation in R_tot. Modeling for a reduced sample size gave broader confidence intervals surrounding the R_tot, R_h and R_a estimates (Figure 10). Estimates of plantation R_tot, R_h and R_a were higher in the random sampling designs than in spatially stratified sampling designs. R_tot and R_a gave higher plantation mean estimates and R_h gave lower plantation mean estimates when rhizosphere data were included in plantation means. Plantation R_tot, and R_a were lower and plantation R_h was higher when samples were taken from the harvest path only (Figure 10; Supplementary Tables S1, S2).

Different spatial patterns were seen in the oil palm plantation for R_tot, R_h, R_a and root density. R_h varied significantly between the different subplots, but R_tot and R_a did not. R_tot and R_a showed significant spatial variation based on the experimental design within the subplots but R_h did not. Manning et al. (2019) took monthly repeated measurements from the different microforms at this site over a year and also found no significant difference in R_h between the different microforms, with variation in water table driving R_h dynamics in the plantation. A neighboring plantation did show significant variation in rates of R_h with surface management microform – this second plantation had an open canopy, unlike the study site here, and the frond piles provided shade that reduced the rates of R_h.

First and foremost, root density was highest next to the palm and decreased as distance from the palm increased. R_tot followed the same pattern as root density, with the highest fluxes measured next to the palm. R_tot fluxes also decreased as distance from palm increased. Similar trends have been seen clearly in other studies (Farmer, 2013; Dariah et al., 2014; Matysek et al., 2018). In this study, there was no significant difference in R_tot with increasing distance from palm after 0.75 m. This was slightly different to root density, which had no significant difference with increasing distance from the palm after 1 m distance from the palm. This study therefore defined the rhizosphere in this plantation as ≤1 m distance from the palm. Modeled estimates of R_a followed the same pattern seen by root density. Modeled estimates of R_h showed no significant variation with distance from palm, remaining relatively constant.

R_tot, R_a and root density varied significantly between the different microforms, suggesting that the spatial patterns in respiration fluxes and root density were at least partially determined by the microforms themselves. There was a clear divide in R_tot, R_a and root density depending on whether the measurements were taken nearer the drainage ditch or the cover plants, with consistently higher measurements nearer the drainage ditches than nearer the cover plants. These trends were also seen in the corresponding frond piles, with higher R_tot, R_a and root density in the frond pile next to the drainage ditch (frond pile-D) and lower R_tot, R_a and root density in the frond pile next to the cover plants (frond pile-C). Furthermore, measurements of R_tot, R_a and root density varied next to the palm in the same transects, showing that the change in root density around the palm began at the palm base.

We propose that the uneven distribution of oil palm roots in space was affected by competition between oil palm and cover plant roots. Competition between oil palm roots and other plants has been shown in a greenhouse experiment in Indonesia (Rahmadhani et al., 2020). Here, oil palm saplings were grown in polybags with herbaceous plants (one plant and one sapling per bag) and root growth was inhibited compared to the polybags that had only an oil palm sapling. Oil palm roots themselves compete for space in mineral soil, when palms are planted 8 m apart (Jourdan and Rey, 1997). In this study, there was minimal root density halfway between two palms, suggesting that roots do not grow as far from the palm in peat soil as they do in mineral soil.

R_tot showed significant spatial variation driven by environmental drivers, as well as driven by distance from palm and surface management microform, which may better explain variation in R_h. Soil pH, soil temperature, air temperature and soil moisture all significantly explained variation in R_tot. Soil pH significantly explained variation in R_tot, with different relationships seen in the different surface microforms. Soil pH varied within the plantation (from 3.08 to 4.45), with pH in soil beneath cover plants being higher than elsewhere. The cover crops consisted of leguminous cover crops, such as Mucuna bracteate. These are planted in oil palm plantations on peat for nitrogen fixation and to preserve soil moisture, in order to minimize the risk of peat subsidence and fires (Othman et al., 2012). In this study, total nitrogen was higher in the cover plant microform than in the other surface microforms, suggesting that some of the nitrogen produced by the cover plants entered the peat soil system and increased the pH locally. Nitrogen fertilization with urea has been shown to increase rates of R_h (Comeau et al., 2016). The increase in R_tot in the cover plants may therefore be caused by the increase in ammonium from nitrogen fixation in the legumes increasing rates of R_h, with soil pH acting as an indicator of the process.

R_tot increased significantly as soil temperature increased. Temperature has been shown to increase the rate of respiration due to an increase in activation energy for biochemical reactions (Lloyd and Taylor, 1994). This can be used as a management strategy - shading tropical peat by 90% has been shown to reduce rates of R_h by 30% (Jauhiainen et al., 2014). The relationship between R_tot and soil temperature did not vary with microform in this study, despite the frond pile and cover plants offering shade. Manning et al. (2019) found that soil temperature and R_h varied between soil management microforms when measured over a year. In longer-term studies at other sites, soil temperature has consistently had a significant effect increasing rates of respiration from peat soil, both when roots were present (Farmer, 2013; Sakata et al., 2015) and absent (Comeau, 2016; Hergoualc’h et al., 2017; Ishikura et al., 2018).

R_tot increased significantly as soil moisture decreased. Soil moisture inhibits R_h by preventing heterotrophic micro-organisms from decomposing the peat, due to the absence of oxygen (Hirano et al., 2012; Mishra et al., 2014; Tonks et al., 2017). Longitudinal data from this site showed that all of the surface microforms had higher R_h rates at lower soil moisture levels, with the frond piles having the strongest effects (Manning et al., 2019). Similar trends have been seen at other oil palm plantations on peat soil (Hergoualc’h et al., 2017; Ishikura et al., 2018; Matysek et al., 2018).

Bulk density was significantly higher in the cover plants than in the harvest path or frond piles. Bulk density was not included in the environmental linear mixed effect model in this plantation due to the bulk density being taken, necessarily (because it disturbs the soils), from slightly different locations to the respiration and root measurements. Melling et al. (2013) found higher R_tot measurements when bulk densities were higher. In this plantation, soil bulk density was significantly higher in the cover plants than from the harvest path or frond piles due to compaction of the harvest path by machinery (Melling et al., 2009). It would be expected that there would therefore be an increase in rates of R_tot and R_h in the cover plants at this plantation due to the increase in bulk density.

The plantation-scale estimates calculated using area-weighted upscaling were decided to give the best estimates of plantation R_tot, R_h and R_a. The best estimate of plantation R_tot in this study was 0.158 ± 0.016 g CO₂-C m⁻² h⁻¹, the best estimates of R_h ranged from 0.112 ± 0.016 to 0.147 ± 0.020 g CO₂-C m⁻² h⁻¹, and the best estimates of R_a ranged from 0.011 ± 0.004 to 0.046 ± 0.026 g CO₂-C m⁻² h⁻¹. These estimates take into account the spatial variation within the oil palm plantation, including scaling to include the proportional area of the drainage ditches. None of these values should be used as annual estimates of R_tot, R_a and R_h fluxes because they were not taken over the year and R_tot and R_a have been shown to have significant temporal variation (Manning et al., 2019).

Two methods were used to partition R_tot into R_h and R_a in order to get these estimates. The lower result comes from the linear regression method, which assumes that the only variation in R_tot is due to variation in R_h and that R_a is fixed. The higher result comes from the distance from palm method, which assumes that the only variation in R_tot is due to variation in R_a and that R_h is fixed. This study has shown that R_tot, root density, R_h and R_a all vary spatially within the plantation. Therefore, the best estimate will lie between these values, with these values providing the upper and lower boundary. We propose that the real value of R_h falls between these two values, i.e., 0.130 ± 0.046 g CO₂-C m⁻² h⁻¹. The best estimate of R_a therefore lies between the estimates of 0.011 ± 0.026 and 0.046 ± 0.004 g CO₂-C m⁻² h⁻¹, i.e., 0.029 ± 0.030 g CO₂-C m⁻² h⁻¹.

Estimates of R_tot, R_h and R_a are similar to other reported estimates of these fluxes in the literature. Published results from chamber measurements of R_tot have a mean of 0.207 ± 0.016 g CO₂-C m⁻² h⁻¹ and range between 0.085 and 0.365 g CO₂-C m⁻² h⁻¹ (Murayama and Bakar, 1996; Melling et al., 2005; Farmer, 2013; Melling et al., 2013; Dariah et al., 2014; Husnain et al., 2014; Sakata et al., 2015; Comeau, 2016; Comeau et al., 2016; Hergoualc’h et al., 2017; Ishikura et al., 2018; Matysek et al., 2018; Manning et al., 2019; Cooper et al., 2020).

Published results from chamber measurements of R_h have a mean of 0.152 ± 0.014 g CO₂-C m⁻² h⁻¹ and range from 0.047 to 0.307 g CO₂-C m⁻² h⁻¹ and published results from chamber measurements of R_a have a mean of 0.088 ± 0.018 g CO₂-C m⁻² h⁻¹ and range from 0.001 to 0.290 g CO₂-C m⁻² h⁻¹ (Farmer, 2013; Melling et al., 2013; Dariah et al., 2014; Husnain et al., 2014; Comeau, 2016; Comeau et al., 2016; Hergoualc’h et al., 2017; Ishikura et al., 2018; Matysek et al., 2018; Manning et al., 2019; Cooper et al., 2020).

Obtaining accurate estimates of R_tot, R_h and R_a is essential in order to create precise flux estimates for climate modeling and for local to global policy decision-making. Published estimates of R_tot, R_h and R_a vary by factors of 4.3, 6.5 and 290.0, respectively. Five reasons why this range is so large include: (1) different partitioning methodologies, (2) different sample sizes, (3) within plantation spatial variations, (4) within plantation micro-climates and (5) variation in temporal dynamics and seasonality. Here we address some of the errors that can be brought in by not taking the first four reasons into account. The fifth reason is explored in more detail in Manning et al. (2019), where significant variation between R_tot and R_a was driven by temporal changes in environmental drivers at this site.

Firstly, this study considered how two different partitioning methodologies gave different results. Estimates of R_h and R_a measured in this study varied by 30 and 61% respectively, considering the straight mean averaging results, and between 24 and 318% for the area-weighted upscaling results. This is important because R_h is often compared between studies without different partitioning methodologies being considered. We recommend using multiple methods to partition R_tot into R_h and R_a where possible, to reduce experimental bias being interpreted as between site variations. We also recommend that modeling studies take partitioning method into account as a covariate.

Secondly, sample size was explored in this research. Sample sizes for R_tot in the literature range from 3 to 72, with an average of 24 samples (Farmer, 2013; Melling et al., 2013; Dariah et al., 2014; Husnain et al., 2014; Comeau, 2016; Comeau et al., 2016; Hergoualc’h et al., 2017; Ishikura et al., 2018; Matysek et al., 2018; Manning et al., 2019; Cooper et al., 2020). This study showed that changing the sample size of the dataset gave different results for R_tot, R_h and R_a. Using resampling techniques to model the sample size of random sampling from 5 to 10, 20 and 35 samples reduced the confidence intervals of the estimate of R_tot, R_h and R_a, making sure that it was more accurate. A power analysis on the dataset in this study suggested that 35 samples were sufficient to give a precise estimate of plantation-scale R_tot, when the samples were stratified based on distance from palm and surface management microform. This suggests that future sampling designs could have more accuracy with larger sample sizes than the average found in the literature.

Thirdly, within-plantation spatial variation was shown to give significant variation in R_tot, R_h and R_a in this study. Therefore, estimates of R_tot, R_h and R_a could be inaccurate if spatial variation is not taken into account. The high root density around the palm and high R_a fluxes in the rhizosphere had a large influence on plantation-scale estimates of respiration. Modeling the difference in R_tot, R_h and R_a in this study when the rhizosphere was excluded, gave 22, 11 and 50% lower estimates of R_tot, R_h and R_a, respectively. Furthermore, estimates of plantation-scale R_tot reduced by 34% when results were scaled up using area-weighted upscaling as opposed to straight mean averaging. This was due to the large contribution from R_a next to the palm that was overrepresented without weighting. Estimates of R_a reduced by 81 and 64% between straight mean averaging and area-weighted upscaling for the results calculated using the distance from palm method and linear regression method, respectively. Similarly, estimates of R_h reduced by 9 and 2% between straight mean averaging and area-weighted upscaling for the results calculated using the distance from palm method and linear regression method. Taken collectively, these results highlight the importance of the spatial variation caused by the rhizosphere in plantation-scale estimates, particularly for accurate estimates of R_tot and R_a.

Fourthly, within plantation micro-climates (aka the surface management microforms) were investigated in this study and one of the key results from this paper was that R_tot and R_a varied between the different surface management microforms outside of the rhizosphere. Here and in Manning et al. (2019) we show different dynamics in soil organic matter mineralization in the different surface management microforms, highlighting the importance of representing these microforms in the plantation-scale respiration estimates. Our sensitivity analysis shows a reduction in plantation-scale R_tot and R_a and an increase in R_h if samples were taken from the harvest path transect only. Other studies measuring R_tot, R_h and R_a from oil palm plantations have focused on measurements in the harvest path, with the exception of Manning et al. (2019), which may lead to an overestimation in plantation R_h.