3.1. Carbon accounting

GHG protocol, initiated by the World Resources Institute and the World Business Council for Sustainable Development, serves as a prominent framework for corporate accounting and reporting emissions. The GHG Protocol Corporate Accounting and Reporting Standard (Protocol and Initiative, Reference Protocol and Initiative2004) categorizes greenhouse gas emissions into three scopes: Scope 1, encompassing direct emissions from owned or controlled sources; Scope 2, involving indirect emissions from purchased energy; and Scope 3, addressing indirect emissions along the entire value chain. We have implemented Application Programming Interface (APIs) for the scalable computation of Scopes 1, 2, and 3 emissions as per the GHG protocol.

3.2. Carbon hotspot detection

The carbon hotspot detection module ingests the carbon footprints of assets along with other relevant associated derived features and selects the prediction and outlier detection algorithm from the library of options, such as linear regressor, decision tree, random forest, Multilayer Perceptron (MLP) regressor, gradient boosting, XGBoost, PyOD, isolation forest, etc., to identify the list of low carbon performing assets/operations. It mainly consists of two sub-modules, namely, (i) emission prediction model and (ii) outlier (note that we use the terms, outlier, anomalous, and hotspot interchangeably) detection model as follows:

3.2.2. Outlier detection model

In this module, we support out-of-the-box outlier detection modules (e.g., PyOD, Isolation Forest, etc.) as well as advanced methods such as multi-dimensional subset scan (MDSS). We used an extension of the MDSS (Neill et al., Reference Neill, McFowland and Zheng2013; Zhang and Neill, Reference Zhang and Neill2016) algorithm for subset scanning of the data. Figure 2 shows the algorithm provided in Zhang and Neill (Reference Zhang and Neill2016). This methodology is able to identify the most anomalous subgroup of feature space in linear time, amongst the exponentially many possible ones, enabling tractable subgroup analysis. The general form of the method is

(3.1) $$ {S}^{\ast }=\mathrm{MDSS}\left(\unicode{x1D53B},\unicode{x1D53C},\mathrm{scor}{\mathrm{e}}_{\mathrm{bias}}\right){S}^{\ast }= MDSS\left(\unicode{x1D53B},\unicode{x1D53C}, scor{e}_{bias}\right) $$

where $ {S}^{\ast } $ is the most anomalous subgroup, $ \unicode{x1D53B}\unicode{x1D53B} $ is a dataset with outcomes $ Y $ and features $ \unicode{x1D54F} $ , $ \unicode{x1D53C} $ are a set of expectations for $ Y $ , and $ \mathrm{scor}{\mathrm{e}}_{\mathrm{bias}} scor{e}_{bias} $ is an expectation-based scoring statistic that measures the amount of anomalousness between subgroup observations and their expectations. The goal of MDSS is to identify a subset of the data $ d(S)\subseteq \unicode{x1D53B} $ corresponding to subgroup $ S $ that maximizes $ \mathrm{scor}{\mathrm{e}}_{\mathrm{bias}} scor{e}_{bias} $ . For $ scor{e}_{bias} $ , we use log-likelihood ratio defined as $ F(S)=\mathit{\log}\left(\mathit{\Pr}\left(D|{H}_1(S)\right)/\mathit{\Pr}\left(D|{H}_0\right)\right) $ . The alternate hypothesis $ {H}_1(S) $ assumes that datapoints $ {x}_i\in d(S) $ are drawn with mean $ q{\mu}_i $ and datapoints $ {x}_i\hskip0.2em \notin \hskip0.2em d(S) $ are drawn from mean $ {\mu}_i $ , for constant multiplicative factor q > 1 known as relative risk. The null hypothesis $ {H}_0 $ assumes that all datapoints, including $ {x}_i\in d(S) $ are drawn with mean $ {\mu}_i $ . This definition satisfies the additive linear-time subset scanning property, which is required for MDSS to be tractable (Speakman et al., Reference Speakman, Somanchi, McFowland and Neill2016). Equation (3.2) gives the general $ \mathrm{scor}{\mathrm{e}}_{\mathrm{bias}}\mathrm{scor}{\mathrm{e}}_{\mathrm{bias}} $ used by MDSS.

(3.2) $$ {\mathrm{score}}_{\mathrm{bias}}(S)=\underset{q>1}{\max}\sum \limits_{x_i\in d(S)}\left(\log \mathit{\Pr}\left({x}_i|q{\mu}_i\right)-\log \mathit{\Pr}\left({x}_i|{\mu}_i\right)\right) $$

Figure 2. MDSS pseudocode.

3.3. Explainability

The interpretation of machine learning models is crucial for understanding feature attribution, especially as many models are inherently not transparent. In this module, we employ the widely recognized model-agnostic methodology known as SHapely Additive exPlanations (SHAP) (Lundberg and Lee, Reference Lundberg, Lee, Guyon, Luxburg, Bengio, Wallach, Fergus, Vishwanathan and Garnett2017) to gain insights into both local and global interpretability for carbon hotspots. SHAP is a game theoretic approach to explain the output of any machine learning model. It connects optimal credit allocation with local explanations using the classic Shapley values from game theory and their related extensions. The SHAP values calculated in our framework represent the marginal contribution of each associated feature to the carbon footprint of an asset. To compute these SHAP values, we leverage the SHAP library available at https://pypi.org/project/shap/. This approach ensures a robust and mathematically sound representation of feature importance. Moreover, our module is designed with a pluggable framework, allowing seamless integration of other explainable models into the workflow.

3.4. Counterfactual queries and recommendation engine

In this module, we provide a framework to leverage the state-of-the-art as well as custom counterfactual queries and recommendation engines. This framework has an inbuilt recommendation module that leverages the diverse counterfactual explanations (DiCE) methodology from Mothilal et al. (Reference Mothilal, Sharma and Tan2020). DiCE generates counterfactual explanations for any ML model through perturbations within a feasible range that change the output of a machine learning model. It also supports simple constraints on features to ensure the feasibility of the generated counterfactual examples. This framework takes input from the user on the set of controllable features with their feasible range and the expected target emissions for intervenable actions and generates the set of best recommendations.

4.1. Overview

Figure 3 gives an overview of palm cultivation. The first step is to produce seedlings from pre-germinated seeds in nurseries, along with land preparation. Then the seedlings are planted manually and fertilizers containing elements such as nitrogen (N), potassium (K), and phosphorus (P) are applied. Fertilizers used are manufactured in factories from where they are taken to the port nearest to the plantations via shipping tankers. Then they are delivered by heavy goods vehicles (HGV) such as trucks to the plantations. Along with fertilizers, diesel is also transported to the plantations, where it is used in irrigation, transportation of fertilizers and workers to the farms, land preparation, and maintenance. It takes about 3–4 years for oil palms to produce fruits suitable for harvest. Palm trees continue to produce fruit for around 30 years and they are harvested periodically. Palm is harvested manually wherein fronds are cut off to dislodge fresh fruit bunches (FFB), which fall to the ground and are then collected. The fruit bunches are transported to oil mills by light good vehicles, where they are subjected to industrial processing to obtain crude palm oil (CPO) and crude palm kernel oil, as well as by-products such as empty fruit bunches (EFB), fibers, shells and palm oil mill effluent (POME). These by-products are returned to the field as manure. Electricity is used in the oil extraction phase to power machinery such as threshers and motors of conveyors.

Figure 3. Palm cultivation: overview.

GHG emissions from the palm supply chain can be divided into the following stages—manufacturing, agriculture, and transportation and electricity usage. Manufacturing includes all the emissions resulting from the production of fertilizers and electricity and the extraction of crude palm oil. Agriculture includes all the emissions resulting from activities related to the planting and harvesting of palm produce. Transportation includes emissions resulting from transporting fertilizers and diesel to palm plantations, and FFB produce from plantations to mills. Lastly, emissions due to electricity usage include all the indirect emissions resulting from the consumption of electricity in the supply chain.

4.2. Palm plantation data

Palm plantations are divided into smaller units called blocks. We use our framework to analyze the performance of palm plantations at the block level, using data from 25 palm blocks for 14 years, from the year of the plantation, up to the 14th year of cultivation of each block. We consider important factors which impact palm cultivation and capture the carbon emission performance of blocks, such as nitrogen content in the fertilizer application, carbon and nitrogen content in manure application (EFBh and pruned fronds), age of the farm block, initial soil organic carbon content at the depth of 0–15 cm and 15–45 cm, annual yield and the weather parameters (annual precipitation and temperature statistics). Data is obtained from a leading producer of palm oil. For the sake of anonymity, we will not disclose the name of the palm oil producer. Along with the farm plantation data, we also use relevant data from the cradle-to-gate extent of the palm oil supply chain for other stages - manufacturing, transportation, and electricity.

Figure 4 shows the heatmap of the annual yield of the 25 blocks at different ages, where the yield of each block has been normalized using min-max scaling. Since palm plantations produce their first harvest in the third–fourth year after plantation, the yield is near zero for the first 3 years across all the blocks. We see that the farm blocks have more yield as their age increases, with maximum yield observed around the 12th year of plantation across most of the blocks.

Figure 4. Heatmap showing normalized annual yield from the farm blocks at different ages.

5.1. Carbon footprint of palm blocks

We have used the physics-based denitrification decomposition (DNDC) model (Gilhespy et al., Reference Gilhespy, Anthony, Cardenas, Chadwick, del Prado, Li, Misselbrook, Rees, Salas, Sanz-Cobena, Smith, Tilston, Topp, Vetter and Yeluripati2014) to estimate the carbon footprint of palm plantations. The model bridges the chemical reactions in soil with the ecological drivers (climate, soil, vegetation, and farming practices) and environmental factors (temperature, moisture, pH, redox potential (Eh), and simulates carbon and nitrogen dynamics. The DNDC model predicts crop growth, soil temperature and moisture, soil carbon sequestration, and emission of CO₂, CH₄, and N₂O along with other trace gases.

Figure 5 shows the inputs, outputs, and components of the model. The inputs to the model are the soil data (type, clay fraction, amount of initial soil organic carbon), weather data (temperature, precipitation, wind speed, solar irradiation), crop data, farming practices (fertilizer, manure, and irrigation schedules, tillage information), planting and harvesting schedules. The two-component DNDC model has six submodules, namely, soil climate, crop growth, denitrification and decomposition, nitrification, and fermentation.

Figure 5. Denitrification-decomposition (DNDC) model workflow (Gilhespy et al., Reference Gilhespy, Anthony, Cardenas, Chadwick, del Prado, Li, Misselbrook, Rees, Salas, Sanz-Cobena, Smith, Tilston, Topp, Vetter and Yeluripati2014).

The soil climate sub-model computes environmental factors within the soil, such as soil temperature, moisture profiles (0–50 cm), pH, and redox potential. These calculations are based on soil texture properties, air temperature, precipitation, and plant water uptake. In parallel, the plant growth sub-model predicts daily water and nitrogen uptake by plants, as well as daily growth, root respiration, and biomass partitioning into grain, stems, and roots. The decomposition model integrates inputs from the soil climate and plant growth models to determine the quantities of substrates, including ammonium ( $ \mathrm{N}{\mathrm{H}}_4^{+} $ ), nitrates ( $ \mathrm{N}{\mathrm{O}}_3^{-} $ ), and dissolved organic carbon.

Depending on the soil moisture profile and redox potential levels, one of the nitrification, denitrification, or fermentation sub-models is activated. Based on the level of soil moisture profile and redox potential, any one of the sub-models of nitrification, denitrification, or fermentation will be invoked. Under aerobic conditions, the nitrification sub-model calculates the nitrification process, leading to the production of $ {\mathrm{N}}_2\mathrm{O} $ , nitric oxide, and ( $ \mathrm{N}{\mathrm{H}}_3 $ ). Conversely, during anaerobic conditions, the denitrification and/or fermentation sub-models estimate the emission of $ {\mathrm{N}}_2\mathrm{O} $ and $ \mathrm{C}{\mathrm{H}}_4 $ , respectively.

The residual crop material integrated into the soil will be categorized into three distinct soil litter pools—namely, highly labile, labile, and resistant litter pools based on their carbon-to-nitrogen (C/N) ratio. This incorporation of litter serves as vital input for the storage of soil organic matter (SOM), thereby connecting the plant and soil components within a biogeochemical system. Within the DNDC framework, SOM is distributed across four primary pools: plant residue (i.e., litter), microbial biomass, humads (i.e., active humus), and passive humus. Each pool comprises two or three sub-pools characterized by varying specific decomposition rates. The daily decomposition rate for each sub-pool is governed by factors such as pool size, specific decomposition rate, soil clay content, nitrogen availability, soil temperature, and soil moisture. As the soil organic carbon (SOC) within a pool undergoes decomposition, a portion of the carbon is released as CO2, while the remainder is allocated to other SOC pools. It is important to note that our study focuses solely on greenhouse gas (GHG) emissions from farming activities, and we do not consider changes in soil organic matter.

This DNDC model has been widely accepted and analyzed by the researchers in last two decades, and several validations with field measurements have been studied for 60+ crops in various countries (Gilhespy et al., Reference Gilhespy, Anthony, Cardenas, Chadwick, del Prado, Li, Misselbrook, Rees, Salas, Sanz-Cobena, Smith, Tilston, Topp, Vetter and Yeluripati2014). For example, a field validation experiment of DNDC model for CH4 and NOx emissions from paddy fields is studied in India and it is reported that the DNDC model satisfactorily simulated the seasonal variations in GHG emissions for different land management (Babu et al., Reference Babu, Li, Frolking, Nayak and Adhya2006).

Our framework computes carbon footprint at each stage from cradle to gate, that is raw material extraction to CPO production leveraging the APIs mentioned in Section 3.1 and the DNDC model. Emissions related to activities within the organizational boundary were computed using Scope1, Scope2 APIs, and DNDC model. Upstream emissions (i.e., raw material extraction and manufacturing of fertilizer) and emissions related to the transportation of upstream raw material and fertilizer were computed using the Scope3 logistic API. At the farming stage, the DNDC tool has been used to precisely compute the carbon footprint due to farming activities. The framework uses “CO₂e” as the standard unit for measuring carbon footprint. For any quantity and type of greenhouse gas, CO₂e signifies the amount of CO₂ which would have the equivalent global warming impact. A quantity of GHG can be expressed as CO₂e by multiplying the amount of the GHG by its global warming potential (GWP), where GWP is the heat absorbed by any greenhouse gas in the atmosphere, as a multiple of the heat that would be absorbed by the same mass of CO₂. For example, the GWP value for N₂0 is provided as 265 in the IPCC Fifth Assessment Report, 2014 (AR5) (Myhre et al., Reference Myhre, Shindell, Bréon, Collins, Fuglestvedt, Huang, Koch, Lamarque, Lee, Mendoza, Nakajima, Robock, Stephens, Takemura and Zhang2013; Myhre et al., Reference Myhre, Shindell, Bréon, Collins, Fuglestvedt, Huang, Koch, Lamarque, Lee, Mendoza, Tignor, Allen, Boschung, Nauels, Xia, Bex and Midgley2014). Therefore, if 1 kg of N₂O is emitted, this can be expressed as 265 kg of CO₂e (1 kg N₂O * 265 = 265 kg CO₂e).

The emissions were divided into four stages—farming, manufacturing, transportation, and electricity as mentioned in Section 4.1. Figure 6 shows the stage-wise distribution of carbon emissions in the cradle-to-gate palm oil supply chain, averaged across the 14 years of plantation data. The plot shows that the farming stage is the major source of emissions in the supply chain. Therefore, in this paper, we will focus on the farming stage, and use our framework for block-level analysis of palm plantation.

Figure 6. Stage-wise distribution of carbon emission.

The DNDC model gives the emission of GHG gases such as CO₂, CH₄, and N₂O from the farming stage at daily intervals. For our analysis, we have aggregated the GHG emissions to annual temporal resolution. Since all the palm plantation blocks are located in the upland mineral soil, the methane emission from the farming stage is negligible. We have considered only the soil CO₂ emission from the decomposition process and N₂O emission from the nitrification and denitrification process to compute the total carbon footprints (CO₂e).

Figure 7 shows the heatmap of annual CO₂e emission of the 25 palm plantation blocks at different ages, where the emission numbers are normalized across all the blocks using min–max scaling. All 25 blocks receive an equal amount of fertilizers for the first 4 years and for subsequent years, the amount of fertilizers and manures are determined based on the age and soil testing. We can see that the annual CO₂e emissions are low and almost in a similar range for the first 4 years and start to show higher value as the age increases across all the blocks.

Figure 7. Normalized annual equivalent carbon emission (CO₂e) from farm blocks at different ages.

5.2. Block-level hotspot identification

The annual carbon footprint of the farm blocks along with the plantation data detailed in Section 4.2 is used for hotspot identification. However, we do not use meteorological variables such as temperature and solar irradiation in identifying the hotspot blocks, even though these variables were used in the DNDC-based estimation of carbon footprint. This is because since the plantation blocks are in close proximity, these variables have similar values across the blocks and, therefore, fail to provide significant discriminatory insights.

To begin, we utilize the emission prediction module to assess the carbon footprint of palm blocks at the block level. We apply a range of well-established machine learning models, such as linear regression, decision tree, random forest, gradient boosting, and MLP regressor. Notably, the random forest regressor exhibits superior performance compared to the other models. The process of determining the optimal hyperparameters for the random forest models is conducted through thorough grid search experiments. This involves methodically examining different hyperparameter values to identify the setup that maximizes performance, ensuring a well-informed and precise model selection.

The predicted carbon emissions serve as the expected emissions $ \unicode{x1D53C} $ , to be used by the MDSS-based outlier detection module. We use the Gaussian scoring function (Speakman et al., Reference Speakman, Somanchi, McFowland and Neill2016) as $ \mathrm{scor}{\mathrm{e}}_{\mathrm{bias}} $ , which uses Gaussian distribution to model the log-likelihood ratio statistic defined in Eqn. (3.2). This serves as the statistical measure of divergence between subgroup observations and their expectations. Gaussian scoring function for a subgroup $ S $ is given by (Speakman et al., Reference Speakman, Somanchi, McFowland and Neill2016):

(5.1) $$ {\mathrm{score}}_{\mathrm{bias}}(S)=\underset{q>1}{\max}\sum \limits_{x_i\in d(S)}{y}_i\hat{\mu}\frac{\left(q-1\right)}{{\hat{\sigma}}^2}+\sum \limits_{x_i\in d(S)}{\hat{\mu}}^2\left(\frac{1-{q}^2}{2{\hat{\sigma}}^2}\right) $$

where, $ {y}_i $ are the observed values in the subset $ d(S) $ belonging to subgroup $ S $ , $ \hat{\mu} $ , and $ {\hat{\sigma}}^2 $ are the expected mean and variance of the subgroup. We perform MDSS-based outlier detection, as outlined in 3.2.2 to obtain hotspot farm blocks, with $ \mathrm{scor}{\mathrm{e}}_{\mathrm{bias}} $ as defined in Eqn. (5.1). B2 and B3 are identified as hotspot blocks across multiple years. In Figure 8 the red region depicts the feature space of the anomalous subgroup $ {S}^{\ast } $ identified by MDSS. From this, we observe that blocks with very low or high fertilizer application, high manure application, and moderate yield are identified as hotspots. However, this does not provide explainable insights or help with identifying opportunities for reducing emissions. To address this, we will next analyze the hotspot blocks using explainability and what-if scenarios.

Figure 8. Anomalous subgroup identified by MDSS.

5.3. Hotspot explainability analysis

To understand the factors that dictate the variation in carbon footprints from palm plantations, we have used the SHAP value to identify the marginal contribution of all the relevant features. We begin with learning the prediction model with farming practices, soil properties, and weather-related parameters as input features for predicting the soil CO₂ and N₂O emissions independently. The individual emission regressor model will give us more insights into the dominating features that influence the respective emissions. Figures 9 and 10 show the holistic summary of feature importance for CO₂ and N₂O emission regressor models, respectively, using the SHAP value.

Figure 9. Global explainability for CO₂.

Figure 10. Global explainability for N₂O.

From Figure 9, we observe that for soil CO₂ emissions, the most dominating features are plantation age, the amount of nitrogen content in the manure and fertilizers, annual yield, annual precipitation, and the initial soil organic carbon (SOC) at the depth of 15 cm–45 cm and 0–15 cm. Since the CO₂ emission is the product of soil decomposition it is clear that the above factors are influencing the rate of the decomposition process. The pruned palm fronds and the EFBs are the main sources of manure applications and have a high carbon-to-nitrogen ratio. Because of the high carbon content, the manure takes precedence over fertilizer in the feature importance plot. The amount of manure and fertilizers tends to increase with age and thus the increase in annual yield. For SOC at the depth of 0–15 cm and 15 cm–45 cm, we see a modest impact with a higher value leading to an increase in CO₂ while the lower value has a neutral impact. This can be attributed to the capacity of soil to hold more carbon without getting saturated. The higher SOC content in the soil tends to have low holding capacity and thus more free carbon to participate in the decomposition process.

Similarly from Figure 10, we see that for N₂O emission, the dominating features are the same as those of CO₂ emission but the order of precedence and the impact are quite different. The nitrogen content in fertilizer takes precedence followed by annual precipitation and then the manure application. The N₂O emission is the product of the nitrification and denitrification process as shown in Figure 5. It is evident from the DNDC model workflow, that the presence of nitrate and ammonium ions with high precipitation tends to be a conducive environment for N₂O emission.

Since the dominating features of soil CO₂ and N₂O emissions are similar, it will be appropriate to investigate the feature importance of total carbon footprint (CO₂e). This will also capture the interdependency among soil CO₂ and N₂O emissions and help the enterprise focus on one metric for the overall carbon footprint reduction. Figure 11 shows the global interpretation of feature importance for the CO₂e emission regressor model using the Shapley value. The order of feature importance is the same as that of soil CO₂ and it shows that the proportion of soil CO₂ is larger than N₂O. Based on the features SHAP value, we infer that in general, the following factors lead to an increase in CO₂e—(i) high fertilizer and manure application results in higher emission. (ii) older blocks tend to have a higher carbon footprint associated with them. (iii) Soil organic carbon at depths 0–15 cm (SOC_45) and 15–45 cm (SOC_15) has a modest impact, with high soil organic carbon leading to an increase in carbon impact, while low SOC has a neutral impact.

Figure 11. Global explainability for CO₂e emission.

We also perform an explainability analysis of the identified hotspot blocks to identify factors behind the low performance of the blocks. We use local explainability SHAP plots as mentioned in Section 3.3. The plot shows which factors contribute to the increase or decrease in CO₂e of the hotspot blocks when compared to the baseline of average carbon footprint. Figures 12 and 13 provide the plots for two hotspot blocks.

Figure 12. Palm plantation block B2 - local explainability.

From Figure 12, we observe that high SOC_15 and SOC_45 content, high manure application result in an increase in emissions, while low fertilizer application lowers the carbon footprint. The block also has a moderate yield which could be due to low fertilizer application. From Figure 13, we can infer that very high fertilizer and manure applications are responsible for high emissions. The block has a good yield performance; however, it has a low carbon performance when compared to other blocks. The observations are in accordance with the insights drawn from Figure 11 earlier.

Figure 13. Palm plantation block B3 - local explainability.

The two explainability plots share some common factors, we can observe that high fertilizer and high manure application negatively impact the carbon performance of the farm blocks.

5.4. Counterfactual recommendation

The recommendation Engine introduced in Section 3.4 is used to generate counterfactual recommendation plans for the hotspot blocks to reduce their carbon footprint. We focused on finding the set of intervenable features which can be perturbed to reduce the CO₂e emission without compromising on the yield produced by palm plantations. Figures 14 and 15 show the intervention plans generated by the engine for the two hotspot blocks. We observe that for both blocks, the fertilizer and manure application are identified as the intervenable areas. For B2, we generated the counterfactual query to produce a set of three best diverse recommendations to restrict CO2e emission within the range of 10 to 25 tonnes of CO₂e. From Figure 14, we observe that reducing fertilizer application by around 10% and manure application by around 30% can help in reducing carbon emissions by roughly 20%, with negligible impact on the yield. Similarly, for B3, we generated the counterfactual query to produce the diverse set of four best recommendations without impacting the yield. From Figure 15, we see that reducing manure application by 25%–30% can reduce the carbon emission by around 10% with negligible impact on the yield for the B3 palm plantation block.

Figure 14. B2 counterfactual recommendations.

Figure 15. B3 counterfactual recommendations.