Herbicide weed control increases nutrient leaching compared to mechanical weeding in a large-scale oil palm plantation

Nutrient leaching in intensively managed oil palm plantations can diminish soil fertility and water quality. There is a need to reduce this environmental footprint without sacrificing yield. In a large-scale oil palm plantation in Acrisol soil, we quantified nutrient leaching using a full factorial experiment with two fertilization rates (260 kg N, 50 kg P, and 220 kg K ha⁻¹ yr⁻¹ as conventional practice and 136 kg N, 17 kg P, and 187 kg K ha⁻¹ yr⁻¹, equal to harvest export, as reduced management) and two weeding methods (conventional herbicide application and mechanical weeding as reduced management) replicated in four blocks. Over the course of 1 year, we collected monthly soil pore water at 1.5 m depth in three distinct management zones: palm circle, inter-row, and frond-stacked area. Nutrient leaching in the palm circle was low due to low solute concentrations and small drainage fluxes, probably resulting from large plant uptake. In contrast, nitrate and aluminum leaching losses were high in the inter-row due to the high concentrations and large drainage fluxes, possibly resulting from low plant uptake and low pH. In the frond-stacked area, base cation leaching was high, presumably from frond litter decomposition, but N leaching was low. Mechanical weeding reduced leaching losses of base cations compared to the conventional herbicide weeding probably because herbicides decreased ground vegetation and thus reduced soil nutrient retention. Reduced fertilization rates diminished the nitrate leaching losses. Leaching of total nitrogen in the mechanical weeding with reduced fertilization treatment (32±6 kg N ha⁻¹ yr⁻¹) was less than half of the conventional management (74±20 kg N ha⁻¹ yr⁻¹), whereas yields were not affected by these treatments. Our findings suggest that mechanical weeding and reduced fertilization should be included in the program by the Indonesian Ministry of Agriculture for precision farming (e.g., variable rates with plantation age), particularly for large-scale oil palm plantations. We further suggest including mechanical weeding and reduced fertilization in science-based policy recommendations, such as those endorsed by the Roundtable for Sustainable Palm Oil association.

1 Introduction

Agricultural expansion is a major driver of tropical deforestation (Geist and Lambin, 2002) which has global impacts on carbon sequestration (Asner et al., 2010; van Straaten et al., 2015; Veldkamp et al., 2020), greenhouse gas regulation (e.g., Murdiyarso et al., 2010; Meijide et al., 2020; Veldkamp et al., 2020), and biodiversity (e.g., Clough et al., 2016). Oil palm is the dominant tree cash crop that replaces tropical forest in Southeast Asia (Gibbs et al., 2010; Carlson et al., 2013) due to its high yields, low production costs, and rising global demand (Carter et al., 2007; Corley, 2009; Grass et al., 2020). Currently, Indonesia contributes 57 % of the global palm oil production (FAO, 2020), which is projected to further expand in the future, threatening the remaining tropical forests (Pirker et al., 2016; Vijay et al., 2016). Forest-to-oil palm conversion is associated with a decrease in soil fertility because of high nutrient export via harvest, reduced rates of soil-N cycling, and decreases in soil organic carbon (SOC) and nutrient stocks (van Straaten et al., 2015; Allen et al., 2015, 2016). Declines in soil fertility promote the dependency on fertilizer inputs and threaten the long-term productivity (Syers, 1997), which could further stimulate the expansion of oil palm production in new areas. Leaching contributes to the reduction in soil nutrient stocks and negatively affects water quality, potentially leading to the eutrophication of water bodies. High loads of nutrients in water bodies due to agricultural expansion and intensification, common in temperate areas (Carpenter et al., 1998), are increasingly reported in humid tropical regions (Figueiredo et al., 2010; Teklu et al., 2018). Because of the high precipitation rates, leaching losses can be substantial in intensively managed plantations in the tropics, although deeply weathered tropical soils also have the capacity to retain large quantities of N and P (Neill et al., 2013; Jankowski et al., 2018). Indeed, nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) can be adsorbed by the anion exchange capacity in the subsoil of highly weathered acidic soils (Wong et al., 1990), whereas P can be fixed to Fe and Al (hydr)oxides common in heavily weathered tropical soils (Roy et al., 2016). Nevertheless, reductions in stream water quality have been reported in oil palm cultivation in Malaysia (Luke et al., 2017; Tokuchi et al., 2019). This illustrates the importance of quantifying nutrient leaching losses in areas with expansive oil palm plantations such as Jambi, Indonesia, one of the hotspots of forest conversion to oil palm in Indonesia (Drescher et al., 2016).

Nutrient leaching losses in oil palm plantations are calculated from water drainage fluxes and solute concentrations (Kurniawan et al., 2018). Despite their relatively low drainage fluxes (as a consequence of high evapotranspiration; Röll et al., 2019; Tarigan et al., 2020), large-scale oil palm plantations typically have high fertilization rates that may result in high nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) concentrations in the soil water and large nitrate leaching losses (e.g., Wakelin et al., 2011; Cannavo et al., 2013). In the leachate, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is accompanied by cations (normally bases) because of its negative charge (Cusack et al., 2009; Dubos et al., 2017), further impoverishing highly weathered tropical soils that are inherently low in base cations (Allen et al., 2016; Kurniawan et al., 2018). High fertilization rates are typically applied to support the high yields of oil palm plantations; however, well-adjusted fertilization rates, for example, to levels that compensate for nutrient export through harvest may create opportunities to reduce nutrient leaching losses while maintaining high productivity.

Herbicides are commonly used for weed control in large-scale oil palm plantations. Herbicides are applied close to the palm stems to reduce competition by weeds for nutrients and water and in the inter-rows to facilitate access during harvest (Corley and Tinker, 2016). Herbicides do not only eradicate aboveground vegetative parts but also remove roots, slowing weed regeneration. Consequently, the use of herbicides for weed control can exacerbate nutrient leaching losses because the absence of ground vegetation reduces the uptake and thus retention of nutrients from applied fertilizers (Abdalla et al., 2019). In contrast to herbicide application, mechanical weeding does not eradicate the roots and allows for the relatively fast regeneration of ground vegetation, which could take up redistributed nutrients and thus reduce leaching losses.

In oil palm plantations, different management zones can be distinguished which have to be taken into account when investigating nutrient leaching losses. Typically, we can identify three contrasting management zones in oil palm plantations: (1) the palm circle, an area of 2 m radius around the palm's stem where the fertilizers are applied and weeded, (2) the inter-row, which is unfertilized and where weed control is less frequent than in the palm circle, and (3) the frond-stacked area, usually every second inter-row, where the pruned senesced fronds are piled up and no weeding or fertilization is done. In each management zone, the extent of nutrient leaching losses depends on the interplay of water fluxes, root uptake, and soil nutrient concentrations. Root uptake, which is related to root density, is high inside the palm circle and lower in the inter-row (Lamade et al., 1996; Jourdan and Rey, 1997). In the palm circle, fertilizers are applied, but also the uptake of water and nutrients is highest (Nelson et al., 2006). Hence, large leaching losses may only occur shortly following fertilization if high drainage fluxes occur, for example, directly following intensive rain showers (Banabas et al., 2008a). The inter-row has higher water input from precipitation than the palm circle because of the lower interception by the canopy (Banabas et al., 2008b). Here, root density and thus root uptake are low, resulting in large water fluxes. However, nutrient leaching may be low in the inter-row because there is no direct fertilizer application. The frond-stacked area receives nutrients from the decomposition of nutrient-rich fronds (Kotowska et al., 2016). Furthermore, mulching with senesced fronds prevents runoff and promotes water infiltration owing to the high macroporosity, a result of high organic matter and biological activity (Moradi et al., 2015). Low canopy interception and high water infiltration may generate high water drainage fluxes, resulting in intermediate nutrient leaching losses in this management zone.

In this study, we aimed to quantify nutrient leaching losses in our experiment that was established in an intensively managed, large-scale oil palm plantation in order to assess whether lower management intensity (i.e., reduced fertilization rates equal to harvest export and mechanical weeding) can reduce leaching losses without affecting yield. We tested the following hypotheses: (1) leaching losses in the palm circle are larger than in other management zones because of direct fertilizer application, (2) leaching losses under herbicide application are higher than mechanical weeding because of the reduced nutrient retention owing to reduced weed growth, and (3) nutrient leaching fluxes under reduced fertilization rates are lower compared to conventional high rates, but yield is not affected. Our study provides the first systematic quantification of leaching losses, an important environmental footprint of oil palm production, by taking into consideration the different management zones and evaluates the effectiveness of alternative management practices on leaching and yield.

2 Materials and methods

2.1 Study area and experimental design

Our study was conducted in a state-owned oil palm plantation in Jambi province, Indonesia (1^∘43^′8^′′ S, 103^∘23^′53^′′ E, 73 m above sea level). Mean annual air temperature is 26.7±1.0 ^∘C, and mean annual precipitation is 2235±385 mm (1991–2011; data from Sultan Thaha airport, Jambi). During our study period (March 2017–February 2018), the mean daily air temperature was 26.3 ^∘C, and annual precipitation was 2772 mm with a dry period between July and October (precipitation less than 140 mm month⁻¹). The soil is highly weathered loam Acrisol soil (Allen et al., 2015), and nutrient inputs from bulk precipitation in the area, measured in 2013, were 12.9 kg N, 0.4 kg P, and 5.5 kg K ha⁻¹ yr⁻¹ (Kurniawan et al., 2018).

This oil palm plantation was established between 1998 and 2002, and the palms were 16–20 years old during our study period. The plantation is mostly located on flat terrain, and it encompasses 2025 ha with a planting density of approximately 142 palms per hectare spaced 8 m apart. The rows between palms are used alternately for harvesting operations and to pile up senesced fronds, which are regularly pruned to facilitate harvesting of fruits; this frond-stacked area covers approximately 15 % of the plantation area. The palm circle, a 2 m radius from the stem where both fertilizers and herbicides are applied, covers 18 % of the plantation. The remaining 67 % we classified as inter-row, which is not fertilized but weeded twice a year.

In November 2016, a factorial management experiment was established with two fertilization rates and two weeding methods (Darras et al., 2019). For fertilizer treatments, the conventional rates were 260 kg N, 50 kg P, and 220 kg K ha⁻¹ yr⁻¹, whereas the reduced rates were 136 kg N, 17 kg P, and 187 kg K ha⁻¹ yr⁻¹. Reduced fertilization rates were established to compensate for nutrient exports via fruit harvest and were assessed by multiplying the nutrient concentrations measured in the fruit bunches with the annual yield. The fertilizer sources were urea (CH₄N₂O), triple superphosphate (Ca(H₂PO₄)₂•H₂O), and muriate of potash (KCl) in granular forms. Fertilizers were applied following the plantation's standard practices: split in two applications per year (in April and October) and spread in a band at approximately 2 m radius from the palm that was raked before fertilizer application. For both fertilizer treatments, we also applied lime (426 kg dolomite ha⁻¹ yr⁻¹; CaMg(CO₃)₂) and micronutrients (142 kg Micro-Mag ha⁻¹ yr⁻¹ with 0.5 % B₂O₃, 0.5 % CuO, 0.25 % Fe₂O₃, 0.15 % ZnO, 0.1 % MnO, and 18 % MgO), as commonly practiced in large-scale plantations on acidic Acrisol soils (Pahan, 2010). Conventional weed control was done using a herbicide (glyphosate), whereas the alternative method was mechanical weeding using a brush cutter; the cut plant materials were left on the ground. Herbicide was applied following the plantation's standard practice: 1.5 L of glyphosate per hectare per year to the palm circle split four times a year and 0.75 L of glyphosate per hectare per year to the inter-row split two times a year. Mechanical weeding was carried out in the same areas and frequencies as herbicide application. This management experiment was comprised of four replicate blocks each with four plots (50 m × 50 m each) assigned to four treatment combinations: conventional rate–herbicide, conventional rate–mechanical weeding, reduced rate–herbicide, and reduced rate–mechanical weeding.

2.2 Soil water sampling

Over the course of 1 year, we collected monthly soil pore water samples using suction cup lysimeters (P80 ceramic, maximum pore size 1 µm; CeramTec AG, Marktredwitz, Germany). We installed the lysimeters in January 2017, randomly choosing two palms per plot and sampling in the three management zones: (1) within in the palm circle 1 m from the palm stem, (2) in the frond-stacked area about 4 m from the palm stem, and (3) in the inter-row approximately 4 m from the palm stem (Fig. A1 in the Appendix). In total, we installed 96 lysimeters (4 treatments × 4 replicates × 2 subplots × 3 management zones). The lysimeters were inserted into the soil to 1.5 m depth so that the soil pore water was collected well below the rooting depth of 1 m, which is common for oil palm plantations on loam Acrisol soils near our study site (Kurniawan et al., 2018). Starting in March 2017, we sampled soil water by applying a 40 kPa vacuum (Kurniawan et al., 2018; Dechert et al., 2005) to the lysimeters. Water samples were collected in dark glass bottles which were stored in a bucket buried in the field. We consider sufficient the 2-month acclimatization of lysimeters before sampling because soil disturbance was minimized and biochemical processes are rapid in tropical soils. During sampling, we transferred once a week the collected water into plastic bottles which were transported to the field station where they were frozen for storage. Soil water collection continued during 1 month until a volume of 100 mL was collected from each lysimeter or until the end of the month. The frozen water samples were transported by air to the University of Göttingen, Germany, where element concentrations were determined. We measured the concentrations of mineral N (${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$), total dissolved N (TDN), and Cl⁻ using continuous flow injection colorimetry (SEAL Analytical AA3; SEAL Analytical GmbH, Norderstedt, Germany), as described in detail by Kurniawan et al. (2018). We calculated dissolved organic N (DON) as the difference between TDN and mineral N. We measured the concentrations of base cations (Na, K, Ca, Mg), total Al, total Fe, total Mn, total S, and total P using an inductively coupled plasma–atomic emission spectrometer (iCAP 6300; Thermo Fisher Scientific GmbH, Dreieich, Germany).

We determined a partial cation–anion charge balance of the major elements (concentrations greater than 0.03 mg L⁻¹) in soil pore water by converting the concentrations to micromoles per liter (µmol_charge L⁻¹). For this, we assumed S to be in the form of sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) and total Al to have a charge of 3⁺. We calculated the combined contribution of organic acids (RCOO⁻) and bicarbonate (${\mathrm{HCO}}_{\mathrm{3}}^{-}$) as the difference between the measured cations and anions (Kurniawan et al., 2018).

2.3 Modeling water drainage

The water balance was modeled using the water sub-model of the Expert-N software, version 5.0 (Priesack, 2005), which was successfully used in previous research to estimate drainage fluxes from different land uses in Indonesia (Dechert et al., 2005; Kurniawan et al., 2018). The model inputs were climate data (solar radiation, temperature, precipitation, relative humidity, and wind speed) and soil (texture, bulk density, and hydraulic functions) and vegetation characteristics (biomass, leaf area index, and root distribution). The climate data were collected from the climatological station in the plantation (described in detail by Meijide et al., 2017), whereas for the oil palm biomass, we used published data from oil palm plantations near our study site (Kotowska et al., 2015). We measured soil bulk density and porosity in the top 10 cm of each management zone in our study site, whereas for the 10–50 cm depth, these were measured in the inter-row. For soil bulk density and porosity at 50–200 cm depth, as well as for soil texture, soil hydraulic parameters (i.e., water retention curve, saturated hydraulic conductivity, and van Genuchten parameters), and root distribution, we used published data from Allen et al. (2015) and Kurniawan et al. (2018), choosing their studied oil palm plantations closest to our study site. The Expert-N water sub-model calculates daily water drainage based on precipitation, evapotranspiration, canopy interception, runoff, and change in soil water storage. Evapotranspiration is calculated using the Penman–Monteith method (Allen, 1998), applying a plant factor of 1.06 (Meijide et al., 2017) and with plant transpiration based on leaf area index (LAI), plant biomass, and maximum rooting depth. The canopy interception is calculated from the percentage of throughfall and the maximum water storage capacity of the canopy. Runoff is calculated from soil texture and bulk density, which determine the water infiltration rate, and from the slope, which was 5 % (Röll et al., 2019). The vertical water movement is calculated using Richards equation based on soil hydraulic functions (Hillel, 1982).

To model the drainage in the different management zones, we used the measured soil bulk density and porosity in the top 10 cm and adjusted other input parameters to simulate differences in water balance in each management zone (Table A1 in the Appendix). For the palm circle, we set the LAI to 3.65, which is the maximum LAI measured at our site (Fan et al., 2015), to simulate high water uptake in the palm circle (Nelson et al., 2006) and the maximum rooting depth of 1 m which is reported for oil palm plantations near our site (Kurniawan et al., 2018). The percentage throughfall in the palm circle was set at 10 %, and the water storage capacity of the oil palm stem was set at 8.4 mm (Tarigan et al., 2018). For the inter-row, we set the LAI and the maximum rooting depth at half the values of the palm circle (1.8 LAI, 50 cm rooting depth) as roots are shallower between palms (Nelson et al., 2006); throughfall was set at 50 %, and the palm stem's water storage capacity was set at 4.7 mm (based on canopy storage capacity reported by Tarigan et al., 2018). For the frond-stacked area, the LAI was set at 0.75, which is half of the minimum measured in the studied plantation (Darras et al., 2019) because understory vegetation is absent in this zone. Values for interception in the frond-stacked area were set at the same values as the inter-row, whereas the runoff was set at 0 (no overland runoff) because mulching with senesced fronds increases water infiltration and prevents runoff (Tarigan et al., 2016).

For validation of the Expert-N water sub-model outputs, we measured weekly soil water matric potential at 30 and 60 cm depths over the study period and compared the measured values with the modeled matric potential. Matric potential was measured by installing a tensiometer (with a P80 ceramic, maximum pore size 1 µm; CeramTec AG, Marktredwitz, Germany) at each depth in each management zone near two palms in two treatments (i.e., conventional rate–herbicide and reduced rate–mechanical weeding) for a total of 12 tensiometers. We summed the modeled daily drainage at 1.5 m depth to get the monthly drainage fluxes which we then multiplied by the element concentrations in soil water to get the monthly nutrient leaching fluxes.

2.4 Soil biochemical characteristics and nutrient retention efficiency

We measured soil biochemical properties in the same sampling locations (Fig. A1) at the following four depth intervals: 0–5, 5–10, 10–30, and 30–50 cm. In each plot, soil samples from the same management zone were pooled to make one composite sample totaling 192 soil samples (4 treatments × 4 replicates × 3 management zones × 4 depths). The samples were air-dried and sieved (2 mm). We measured pH on a 1 : 4 soil-to-water ratio and effective cation exchange capacity (ECEC) by percolating the soils with 1 mol L⁻¹ of unbuffered NH₄Cl and analyzing the cations (Ca, Mg, K, Na, Al, Fe, Mn) in percolates using ICP-AES (iCAP 6300 Duo view ICP Spectrometer, Thermo Fisher Scientific GmbH, Dreieich, Germany). A subsample was finely ground and analyzed for organic C and total N using a CN analyzer (Vario EL Cube; Elementar Analysis Systems GmbH, Hanau, Germany) and for ¹⁵N natural abundance signature using an isotope ratio mass spectrometer (IRMS; Delta Plus, Finnigan MAT, Bremen, Germany). We calculated the soil element stocks for each depth by multiplying the element concentration with the measured bulk density and adding them for the top 50 cm; other soil characteristics (e.g., pH, ECEC, base saturation) in the top 50 cm of soil were calculated as the depth-weighted average of the sampled depths.

In addition, we calculated the N and base cation retention efficiency in the soil for each experimental treatment and management zone following the following formula: nutrient retention efficiency equals 1 minus the nutrient leaching loss divided by the soil-available nutrient (Kurniawan et al., 2018). We used the gross N mineralization rates in the top 5 cm depth (Table A2) as an index of soil-available N, whereas soil-available base cations were the sum of the stocks of K, Na, Mg, and Ca in the top 10 cm depth (expressed in mol_charge m⁻²).

2.5 Statistical analyses

For soil biochemical properties measured once, we tested for differences among management zones and experimental treatments for the entire 50 cm depth using the analysis of variances (ANOVA) with the Tukey honestly significant difference (HSD) test as a post hoc test. The soil variables that showed non-normal distribution or unequal variances were log-transformed prior to the analysis with the Shapiro–Wilk and Levene's tests, respectively. Base cation and N retention efficiencies were also tested for differences between experimental treatments in the same way. For repeatedly measured variables, i.e., soil pore water solute concentrations and leaching fluxes, we used linear mixed effects (LME) models to assess the differences among management zones and treatments. Management zones and treatment, given by fertilization, weeding, and their interaction, were the fixed effects of the LME, and the random effects were sampling months and replicate plots, which were also nested with management zones and in turn nested with subplots. Differences were assessed with ANOVA followed by the Tukey HSD test (Hothorn et al., 2008). In addition, we included a variance function that allows different variances of the management zones. We also used LME to assess differences in soil water matric potential among management zones with management zone as the fixed effect and measurement day and depth nested with treatment as random effects. Comparability between modeled and measured soil water matric potential for each depth in each management zone (n=50 field measurements) was assessed using the Pearson correlation test. All tests were considered significant at P≤0.05 except for soil pH, for which there was a marginal significance at P=0.06. All statistical analyses were performed with R version 3.6.1 (R Core Team, 2019).

3 Results

3.1 Soil biochemical properties and water balance

Soil biochemical properties in the top 50 cm did not differ between experimental treatments (all P>0.05) but strongly differed among management zones (Table 1). The frond-stacked area, where senesced fronds were regularly piled like mulch material, had higher SOC and total N stocks (P<0.01) compared to the other management zones. The inter-row, with regular weeding but without direct fertilizer and lime inputs, showed lower exchangeable base cation contents (i.e., Ca, Mg, K) compared to the other management zones (P≤0.02) and higher exchangeable Al content than the palm circle (P=0.01). This was reflected in the lower base saturation and higher Al saturation in the inter-row compared to the other zones (P<0.01). Also, inter-row had the lowest ECEC (P<0.01) and marginally lower pH than the palm circle (P=0.06). The palm circle, where fertilizers and lime were applied, had generally comparable exchangeable element contents with the frond-stacked area except for K, which was higher in the palm circle (P<0.01), and for Mn, which was higher in the frond-stacked area (P<0.01).

Table 1 Soil physical and biochemical characteristics (mean ± standard errors, n=4 plots) in the top 50 cm depth for each management zone, averaged across experimental treatments. Means within a row followed by different letters indicate significant differences among management zones (one-way ANOVA with the Tukey HSD or Kruskal–Wallis H test with multiple comparison extensions at P≤0.05). Bulk density measured in the top 10 cm of soil, whereas all the other parameters are for 0–50 cm soil depth: element stocks are the sum of the sampled soil depths (0–5, 5–10, 10–30 and 30–50 cm) and the rest are depth-weighted averages, calculated for each replicate plot. ECEC, effective cation exchange capacity.

There were strong positive correlations between field-measured and modeled soil water matric potential (Fig. 1). The matric potential was generally lowest in the palm circle, intermediate in the inter-row, and highest in the frond-stacked area (P<0.01). This pattern was also reflected in the low drainage flux in the palm circle and high drainage flux in the frond-stacked area (Table 2; Fig. 2). In the palm circle, the low drainage flux resulted from high plant transpiration and interception, whereas the high drainage flux in the frond-stacked area was due to low evapotranspiration and runoff with the senesced frond mulch (Table 2). Compared to annual precipitation, the calculated annual evapotranspiration was 51 %, 31 %, and 38 % in the palm circle, frond-stacked area, and inter-row, respectively; annual drainage fluxes at 1.5 m depth were 20 % of precipitation in the palm circle, 65 % in the frond-stacked area, and 43 % in the inter-row. Over the course of 1 year, the monthly drainage fluxes displayed two peaks, in May and November, which occurred following several days of moderate rainfall. The lowest drainage fluxes were measured during the end of the dry season (Fig. 2).

Table 2 Annual water balance simulated from March 2017 to February 2018 for each management zone.

Figure 1 Pearson correlation test between modeled (red line) and field-measured (black points) soil water matric potential (n=50 field measurements over 1 year) for each management zone at 30 and 60 cm depths.

Figure 2 Monthly water drainage at 1.5 m depth, simulated in each management zone, and daily rainfall from March 2017 to February 2018. The gray shaded areas represent the dry season (precipitation less than 140 mm month⁻¹).

3.2 Differences in leaching losses among management zones and treatments

The frond-stacked area had generally lower ionic charge concentrations compared to the other management zones (Fig. 3). Dominant cations in leachate were Al³⁺, Ca²⁺, Mg²⁺, K⁺, and Na⁺ across experimental treatments and management zones. Dissolved Al concentrations were higher in the inter-row than in the other zones (P<0.01). The Ca²⁺ concentrations were similar in the palm circle and frond-stacked areas (P=0.96), and both were higher than in the inter-row (P<0.01). The concentrations of Mg²⁺ and K⁺ were higher in the palm circle than in the other two management zones (P<0.03). The Na⁺ concentrations were higher in the palm circle and inter-row than in the frond-stacked area (P<0.01). As for N, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations were higher in the inter-row than in the other zones (P<0.01). Across treatments, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was 4 %–18 % of TDN, whereas DON was only 1 %–7 % of TDN. Thus, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was the main form of dissolved N, which was higher in the inter-row (P<0.01). The dominant anion was Cl⁻ with higher concentrations in the palm circle than in the other zones (P<0.01).

Figure 3 Partial cation–anion charge balance of the major solutes (with concentrations greater than 0.03 mg L⁻¹) in soil water at 1.5 m depth for each experimental treatment in the different management zones. The combined concentrations of organic acids (RCOO⁻) and carbonates (${\mathrm{HCO}}_{\mathrm{3}}^{-}$) are calculated as the difference between the measured cations and anions. Treatments: ch is conventional fertilization–herbicide, cw is conventional fertilization–mechanical weeding, rh is reduced fertilization–herbicide, and rw is reduced fertilization–mechanical weeding.

Monthly leaching fluxes showed a common pattern among the major solutes (Fig. 4): two peaks of leaching losses (May and November) following fertilizer applications, whereas lower leaching losses occurred during the dry season from July to October. Leaching fluxes of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and total Al followed a similar spatial pattern as their concentrations, being higher in the inter-row than the other zones (P<0.01; Fig. 4). In contrast, Ca leaching fluxes were higher in the frond-stacked area than the palm circle and inter-row (P<0.01; Fig. 4). The leaching of K, Mg, and Na was comparable among zones (P>0.16; Fig. 4).

Figure 4 Monthly leaching losses at 1.5 m depth (mean ± standard errors, n=4 plots) for each management zone. Black arrows indicate fertilizer applications, and the gray shaded areas represent the dry season (precipitation less than 140 mm month⁻¹).

The reduced intensity of management influenced nutrient leaching losses (Fig. 5; Table 3). Reduced fertilization rates diminished ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching fluxes compared to conventional fertilization (P=0.05; Fig. 5). Mechanical weeding lowered leaching losses of Ca, K, and Mg compared to herbicide weeding (all P≤0.05; Fig. 5).

Table 3 Annual leaching losses at 1.5 m depth for each experimental treatment from March 2017 to February 2018. Values are area-weighted averages of leaching losses in each management zone (mean ± standard errors, n=4 plots). Treatments: ch is conventional fertilization–herbicide, cw is conventional fertilization–mechanical weeding, rh is reduced fertilization–herbicide, and rw is reduced fertilization–mechanical weeding. DON is dissolved organic N, and TDN is total dissolved N.

Figure 5 Average monthly leaching losses at 1.5 m depth for each experimental treatment from March 2017 to February 2018 (means ± standard errors, n=4 plots). Effects of fertilization and weeding were evaluated using linear mixed effects models with monthly sampling as the random effect and plots nested by management zone with subplot as the random effect; there was no significant effect of fertilization and weeding interaction. Treatments: ch is conventional fertilization–herbicide, cw is conventional fertilization–mechanical weeding, rh is reduced fertilization–herbicide, and rw is reduced fertilization–mechanical weeding.

3.3 Annual leaching losses and nutrient retention efficiency

Annual leaching losses in the plantation were calculated from the weighted-average of annual leaching in each management zone: 18 % for the palm circle, 15 % for the frond-stacked area, and 67 % for the inter-row. In proportion to the applied fertilizer, annual leaching losses of TDN (Table 3) were 28 % of the applied N in the herbicide treatment for both conventional and reduced fertilization rates, 24 % in the mechanical weeding with conventional fertilization, and only 19 % in the mechanical weeding with reduced fertilization. The annual leaching of K (Table 3) was 4 % of the applied K fertilizer in the herbicide treatment and 3 % in the mechanical weeding for both fertilization rates. In these highly weathered Acrisol soils with a high capacity for P fixation by Fe and Al (hydr)oxides, we detected no leaching of dissolved P (Table 3).

Both N and base cation retention efficiencies were generally lower in the inter-row compared to the other management zones (P≤0.03), except for reduced fertilization–mechanical weeding for which there were no differences among management zones (Table 4). The area-weighted average N retention efficiency was comparable among experimental treatments (P=0.89), but there was a trend of increasing efficiency with decreasing management intensity (Table 4). Base cation retention efficiency showed strong differences among experimental treatments for each management zone: in the palm circle, it was highest in mechanical weeding and lowest in the herbicide treatment (P=0.04); in the frond-staked area and inter-row, it was lowest in the most intensive management treatment (conventional fertilization–herbicide) and highest in either mechanical weeding or reduced fertilization (P≤0.05; Table 4). The area-weighted averaged base cation retention efficiency was also influenced by weeding method, being both lowest in herbicide treatment and highest in mechanical weeding with conventional fertilization (P=0.03; Table 4).

Table 4 N and base cation retention efficiencies in the soil for each management zone and experimental treatment (means ± standard errors, n=4 plots). Means followed by different lowercase letters indicate differences among experimental treatments for each management zone, whereas different uppercase letters indicate differences among management zones for each experimental treatment (one-way ANOVA with the Tukey HSD or Kruskal–Wallis H test with multiple comparison extensions at P≤0.05). Weighted average is based on the areal coverage of each management zone: 18 % for palm circle, 15 % for frond-stacked, and 67 % for inter-row areas. Treatments: ch is conventional fertilization–herbicide, cw is conventional fertilization–mechanical weeding, rh is reduced fertilization–herbicide, and rw is reduced fertilization–mechanical weeding. See Sect. 2.4 for calculations of N and base cation retention efficiency.

4 Discussion

4.1 Water model and temporal pattern of nutrient leaching losses

To our knowledge, our study is the first that has modeled water drainage fluxes from the different management zones of an oil palm plantation, which makes a comparison with other published values challenging. Modeled annual transpiration rates in the palm circle (Table 2) were remarkably similar to the values estimated with the eddy covariance technique in the same oil palm plantation (827–829 mm yr⁻¹; Meijide et al., 2017; Röll et al., 2019). Furthermore, our average daily transpiration rate (2.3 mm d⁻¹) was within the range of rates measured with drone-based photogrammetry (3±1 mm d⁻¹; Ahongshangbam et al., 2019) performed in the same plantation. The modeled annual runoff in the palm circle and inter-row (Table 2) was also within the range of runoff estimates in oil palm plantations in Jambi province (10 %–20 % of rainfall; Tarigan et al., 2016) and in Papua New Guinea (1.4 %–6 % of rainfall; Banabas et al., 2008b). Considering the areal proportions of the three management zones, the weighted-average drainage flux (1161 mm yr⁻¹) was lower than the estimate for smallholder oil palm plantations near our study site (1614 mm drainage flux with 3418 mm precipitation measured in 2013; Kurniawan et al., 2018). However, higher evapotranspiration rates in large-scale compared to smallholder oil palm plantations in our study area (Röll et al., 2019) may explain these differences. Nevertheless, ratios of drainage flux to annual precipitation were comparable between our study and the study by Kurniawan et al. (2018). We conclude from these comparisons with literature values and on the good agreement between modeled and measured soil water matric potential (Fig. 1) that our modeled water drainage fluxes were reliable. The frond-stacked areas had larger drainage fluxes caused by a combination of low evapotranspiration and runoff (Table 2) and enhanced porosity (indicated by lower bulk density; Table 1) from organic matter that facilitates water infiltration (Moradi et al., 2015). This suggests that piling senesced fronds may amend groundwater recharge which, in turn, could moderate discharge fluctuations in water catchments of oil-palm-converted areas (Tarigan et al., 2020).

The temporal peaks of nutrient leaching fluxes (May and November; Fig. 4) likely resulted from the combined effect of high drainage flux and fertilizer application. Large drainage fluxes might have stimulated the downward transport of nutrients and decreased their residence time in the soil and thus their adsorption into the soil exchange sites (Lohse and Matson, 2005). Large drainage fluxes usually dilute the nutrient concentrations in the soil pore water; however, the combined fertilizer and lime applications were able to maintain high nutrient concentrations as manifested by the parallel peaks of drainage and nutrient leaching fluxes (Figs. 2 and 4). The high ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching following urea N fertilization (Fig. 4) suggests rapid nitrification (Silver et al., 2005), fast ${\mathrm{NO}}_{\mathrm{3}}^{-}$ transport through the soil column, and limited anion adsorption capacity (Wong et al., 1990). The latter was possibly affected by the added Cl⁻ from fertilization with KCl (Fig. 3) which may have saturated the soil anion exchange sites, particularly in this mature plantation which has been intensively fertilized for 16–20 years. Due to its negative charge, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching fluxes are always accompanied by comparable leaching fluxes of positive cations (Dubos et al., 2017; Kurniawan et al., 2018), resulting in similar temporal leaching patterns (Fig. 4). Our findings illustrate that fertilization should be avoided during periods of high drainage fluxes which are related to extended periods of moderate rainfall (Fig. 2). However, it is expected that the reliable prediction of periods with high rainfall and drainage will become even more difficult with climate change, which is increasing uncertainties in rainfall intensity and distribution (Chou et al., 2013; Feng et al., 2013). Fertilization during the dry season is also not advisable because plant uptake is low during this period (Corley and Tinker, 2016) and the application of urea, together with lime, will cause urea to volatilize easily even in these acidic soils (Goh et al., 2003; Pardon et al., 2016).

Our results suggest that there are several viable options to reduce leaching losses without sacrificing production. Spreading fertilizer applications over a longer period and reducing fertilization rates, for example, at a compensatory level equal to harvest export, as we tested in our experiment, are recommendable alternatives to present practices. In addition, the use of organic amendments, such as empty fruit bunches, compost, palm oil mill effluent, or slow-release fertilizers, which have been shown to reduce N leaching in tropical cropping systems (Nyamangara et al., 2003; Steiner et al., 2008; Mohanty et al., 2018), will also reduce leaching losses. Organic fertilizers have the additional advantage of improving soil fertility in oil palm plantations (Comte et al., 2013; Boafo et al., 2020), as was also shown by the mulching of senesced oil palm fronds (i.e., high SOC, total N, ECEC, and base saturation in the frond-stacked area; Table 1).

4.2 Leaching losses in the different management zones

A surprising result, in contrast to our first hypothesis, was that nutrient leaching losses among management zones were generally large in the inter-row, especially for mineral N (largely NO₃; Figs. 3 and 4). We did not expect this because the inter-row did not receive direct fertilizer inputs (see Sect. 2.1). Our results suggest that mineral N was transported via surface and/or subsurface lateral flows from the fertilized palm circle to the inter-row which were only 3 m apart (Fig. A1). We expect that the contribution of surface transport of mineral N was a minor process at our site because of the low runoff (Table 2). In an oil palm plantation in Papua New Guinea as well, the loss of N fertilizer via surface runoff was only 0.3–2.2 kg N ha⁻¹ yr⁻¹ (Banabas et al., 2008b). The dominant form of transport of mineral N in our experiment was likely by subsurface lateral flow. Acrisol soils are characterized by clay translocation to a subsurface soil horizon that can create a stagnating layer above which lateral water flow can occur (Elsenbeer, 2001). Indeed, the clay contents of the Acrisol soils at our study area increase with depth, and soil bulk density at 100–150 cm was greater than at 150–200 cm depth (Allen et al., 2016). In addition, palm roots spreading from the palm circle to the inter-row may create channels for the subsurface lateral flow of dissolved ions such as ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Li and Ghodrati, 1994). Higher mineral N leaching in the inter-row than palm circle had also been observed in a study in Brazil where it was attributed to lower root density and higher N mineralization at increasing distance from the palm's stem (Schroth et al., 2000). Hence, a combination of lower root uptake, higher N mineralization, and subsurface lateral transport (particularly for ${\mathrm{NO}}_{\mathrm{3}}^{-}$) all may have contributed to higher mineral N leaching losses in the inter-row than the palm circle. In the inter-row, the main cation that accompanied the leached ${\mathrm{NO}}_{\mathrm{3}}^{-}$ was Al³⁺ (Figs. 3 and 4). This is because this zone's soil pH (Table 1) was within the Al buffering range (pH 3–5; van Breemen et al., 1983) as this zone had no direct lime application and consequently had a low base saturation (Table 1). Our findings also show that if leaching is measured only within the palm circle, this could lead to a substantial underestimation of mineral N and Al leaching losses.

Despite the direct application of fertilizer, the palm circle had relatively low N leaching losses (Figs. 3 and 4), which was probably due to the large root density facilitating an efficient nutrient uptake (Edy et al., 2020; Nelson et al., 2006). The dominant anion in soil pore water in the palm circle was Cl⁻ (Fig. 3) which was enhanced by the applied KCl fertilizer and which was accompanied by high base cation concentrations relative to dissolved Al (Fig. 3). The former was due to the applied Micro-Mag fertilizer and dolomite (Sect. 2.1) which increased pH and exchangeable bases and rendered Al into an insoluble form (Table 1; Schlesinger and Bernhardt, 2013). Despite their high concentrations, base cation leaching fluxes in the palm circle (Fig. 4) were constrained by the low water drainage flux (Table 2).

Although the frond-stacked area was at the same distance from the palm circle as the inter-row (Fig. A1), mineral N leaching losses were substantially lower (Figs. 3 and 4). Decomposition of nutrient-rich fronds (Kotowska et al., 2016) resulted in high SOC and N stocks (Table 1), which can support a large microbial biomass in this zone (Haron et al., 1998). Immobilization of mineral N by the large microbial biomass converting mobile ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to less mobile organic N may have caused the low mineral N leaching in the frond-stacked area (e.g., Corre et al., 2010). In addition, palm root uptake of nutrients (including mineral N) may have been higher in the frond-stacked area than in the inter-row because roots tend to proliferate in nutrient-rich zones (Table 1; Hodge, 2004). Indeed, studies have shown higher root density and higher water uptake under the frond piles compared to the inter-row (Nelson et al., 2006; Rüegg et al., 2019). The larger base cation leaching in the frond-stacked area compared to the inter-row (Fig. 4) was probably a reflection of the high ECEC, base saturation, and pH in the frond-stacked area (Table 1). These favorable soil characteristics were probably caused by the release of nutrients from the decomposition of frond litter, which contains high base cation concentrations (Kotowska et al., 2016). Finally, the low Al leaching in the frond-stacked area (Figs. 3 and 4) can be explained by the higher soil pH (Table 1). Our results highlight the benefits of piling senesced fronds on the soil to reduce leaching of mineral N and Al, which could otherwise affect ground water quality. In other areas such as Borneo, oil palm plantations were reported to practice the piling of senesced fronds on every inter-row (Rahman et al., 2018). In our study region, this is rarely practiced because it hinders access to palms during harvest. Maybe chopping-up senesced leaves with a shredder before spreading them on the soil can both improve access and at the same time enhance nutrient management of oil palm plantations.

4.3 Leaching losses under a different intensity of management

Management intensity treatments affected nutrient leaching losses with generally lower leaching fluxes under less-intensive management (Fig. 5; Table 3). In line with our second hypothesis, mechanical weeding had lower leaching fluxes of Ca, K, and Mg than the herbicide application (Fig. 5; Table 5). Plots with mechanical weeding had higher ground vegetation cover (Darras et al., 2019) and higher nutrient retention efficiency than herbicide weeding (Table 4). Leaching losses were probably retained better by the faster regrowth of understory vegetation under mechanical weeding. This is in line with studies in temperate forests and in a cedar plantation showing that understory vegetation can take up excess ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in the soil (Olsson and Falkengren-Grerup, 2003) and reduce ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching and the mobilization of Ca and Mg (Fukuzawa et al., 2006; Baba et al., 2011). Denser understory vegetation in oil palm plantations may also positively impact biodiversity by increasing the richness of plant species and soil macrofauna diversity and abundance (Ashton-Butt et al., 2018; Luke et al., 2019), which may facilitate nutrient uptake and recycling. Following the first 3 years after the establishment of the experiment, oil palm yield was approximately 30 Mg of fresh fruit bunches per hectare per year and did not differ among experimental treatments (Fig. A2; Darras et al., 2019). This attests that during the first 3 years, the reduced management intensity did not affect productivity. However, long-term monitoring of yield is essential as it may take a longer period before the yield responds to our experimental treatments (e.g., Tao et al., 2017). The costs of the two weeding treatments (i.e., herbicide vs mechanical) were not different because it is a common practice to combine the use of herbicide with the periodic mechanical cutting of resistant ground vegetation (Pahan, 2010; Darras et al., 2019). In addition, the use of glyphosate has been associated with possible health risks to workers and the environment (van Bruggen et al., 2018). In summary, our results advocate for a more sustainable management with mechanical weeding compared to herbicide application.

The decrease in ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching with reduced N fertilization rates without affecting yield supports our third hypothesis. Our results suggest that excess N applied with the conventional fertilization rate (above harvest export; Sect. 2.1) was largely lost through leaching (Table 3) as there were no differences in total N stocks (Sect. 3.1), mineral N levels (Darras et al., 2019), N retention efficiency (Table 4), and oil palm yield (Darras et al., 2019). The reduced K fertilization did not affect K leaching (Fig. 5) probably because K fertilization rates were only reduced by 15 % of the conventional rate owing to high K export with harvested oil palm fruits (Sect. 2.1). Our study provides evidence that this mature (16–20 years old) plantation with conventional management was over-fertilized with N, and we suggest that the inclusion of lower N fertilization rates (related to N export with fruit bunches) in the Indonesian program for precision farming (Ministry of Agriculture of Indonesia, 2016) will substantially and quickly improve the environmental footprint of oil palm production.

Compared to other fertilized tropical plantations (Table A3), our plantation had similar N leaching estimates to those reported in another oil palm study using a model validated with field data from Sumatra, Indonesia (Pardon et al., 2020). In contrast, lower N leaching losses were reported in other large-scale oil palm plantations in similar soils with comparable fertilization rates (Omoti et al., 1983; Tung et al., 2009). However, in these studies, leaching losses were exclusively measured in the palm circle (Omoti et al., 1983), or the sampling location was not specified (Tung et al., 2009). Both studies may thus have underestimated N leaching because our results showed the highest contribution to leaching losses were from the inter-row (Figs. 3 and 4). N leaching fluxes in our plantation were also higher than fluxes reported from smallholder oil palm plantations in the same area owing to their lower fertilization rates (Kurniawan et al., 2018). In contrast, N leaching in our plantation was lower than from an oil palm plantation or coffee agroforestry system in volcanic soils (Banabas et al., 2008b; Tully et al., 2012; Cannavo et al., 2013). This may be caused by the inherently higher nutrient content and high porosity of these volcanic soils that facilitate high infiltration rates. N leaching losses from our plantation were also lower compared to banana plantations, which had substantially higher fertilization rates (Wakelin et al., 2011; Armour et al., 2013).

The high fluxes of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and Al at 1.5 m depth imply a substantial risk of groundwater pollution. During the period of high drainage fluxes following fertilization, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations in soil pore water reached concentrations of 20–40 mg ${\mathrm{NO}}_{\mathrm{3}}^{-}$ L⁻¹ in the inter-row (covering 67 % of the plantation area), which is close to the upper limit of 50 mg ${\mathrm{NO}}_{\mathrm{3}}^{-}$ L⁻¹ for drinking water (WHO, 2011). Al concentrations in soil pore water even exceeded the limit of 0.2 mg Al L⁻¹ in 60 % of the samples. This does not automatically mean that surface water will be contaminated as ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and Al concentrations can be diluted and partially retained in the soil (Harmand et al., 2010; Jankowski et al., 2018) or denitrified (Wakelin et al., 2011). Such processes are especially effective in riparian buffers, which can mitigate the transport of these agricultural pollutants to streams (Luke et al., 2017; Chellaiah and Yule, 2018). Our results thus support the importance of restoring riparian buffers in areas converted to oil palm plantations, which is also an important sustainability criterion endorsed by the Roundtable for Sustainable Palm Oil (RSPO) association (RSPO, 2018) that may provide additional regulation services (Woodham et al., 2019).

5 Conclusions

Our findings show that nutrient leaching losses in an oil palm plantation differed among management zones as a result of fertilization, liming, mulching, and different drainage fluxes. The implementation of mechanical weeding with reduced fertilization rates was effective in reducing nutrient leaching losses without affecting yield during the first 3 years of this experiment. The long-term investigation of this management experiment is important and has been planned in order to get a reliable response of yield and to make a more holistic economic analysis that includes the valuation of regulation services. Greenhouse gas emissions should also be quantified as another important parameter of the environmental footprint of oil palm production. Our ultimate goal is that our present and future findings will be incorporated into science-based policy recommendations such as those endorsed by the RSPO.

Appendix A

Table A1 Parameters used in the Expert-N water sub-model for each management zone.

Table A2 Gross N mineralization rates (mean ± standard errors, n=4 plots) in the top 5 cm of soil for each treatment and management zone in a large-scale plantation in Jambi, Indonesia. Measurements were done on intact soil cores in February 2018 using the ¹⁵N pool dilution technique, as described in detail by Allen et al. (2015). Treatments: ch is conventional fertilization–herbicide, cw is conventional fertilization–mechanical weeding, rh is reduced fertilization–herbicide, and rw is reduced fertilization–mechanical weeding.

Note that these data are not included in the main paper to avoid redundant publication as they were already included in another paper presently in review.

Table A3 Literature comparison of annual N fertilization and total N leaching losses across tropical plantations.

Figure A1 Lysimeter locations at each treatment plot with two subplots (blue rectangles) that each included the three management zones (blue crosses): (1) lysimeters in the palm circle were 1 m from the palm stem, (2) in the frond-stacked area about 4 m from the palm stem, and (3) in the inter-row approximately 4 m from the palm stem.

Figure A2 Annual yield of each experimental treatment from 2017 to 2019. Treatments: ch is conventional fertilization–herbicide, cw is conventional fertilization–mechanical weeding, rh is reduced fertilization–herbicide, and rw is reduced fertilization–mechanical weeding. Note that yield was measured by weighing the harvested fresh fruit bunches from each palm in the inner 30 m × 30 m area of each plot.