^{1}

^{2}

^{*}

^{3}

^{4}

^{5}

^{1}

^{2}

^{2}

^{6}

^{7}

^{8}

^{9}

^{10}

^{1}

^{1}

^{11}

^{12}

Climate change is being intensified by anthropogenic emission of greenhouse gasses, highlighting the value of forests for carbon dioxide storing carbon in their biomass. Seasonally dry tropical forests are a neglected, threatened, but potentially critical biome for helping mitigate climate change. In South America, knowing the amount and distribution of carbon in Caatinga seasonally dry vegetation is essential to understand its contribution to the global carbon cycle and subsequently design a strategic plan for its conservation. The present study aimed to model and map the spatial distribution of the potential forest biomass stock across 32 forest fragments of Caatinga, in the state of Bahia, Brazil, using regression kriging and Inverse Square of Distance techniques, building from point measurements of vegetation biomass made on-the-ground in ecological plots. First, a model for estimating biomass was fitted as a function of environmental variables to apply regression kriging, and then applied to the maps of the selected components. Elevation, temperature, and precipitation explained 46% of the biomass variations in the Caatinga. The model residuals showed strong spatial dependence and were mapped based on geostatistical criteria, selecting the spherical semivariogram model for interpolation by ordinary kriging. Biomass was also mapped by the Inverse Square of Distance approach. The quality of the regression model suggests that there is good potential for estimating biomass here from environmental variables. The regression kriging showed greater detail in the spatial distribution and revealed a spatial trend of increasing biomass from the north to south of the domain. Additional studies with greater sampling intensity and the use of other explanatory variables are suggested to improve the model, as well as to maximize the technique’s ability to capture the actual biomass behavior in this newly studied seasonally dry ecosystem.

Seasonally dry tropical forests (SDTFs) act as important CO_{2} sinks as they promote carbon fixation from the atmosphere in their plant biomass and contribute to the quantification and cycling of nutrients (

While the importance of this information for the conservation, management, and sustainable use of SDTFs is recognized (

Measuring environmental variables is an essential step towards understanding on how to better quantify ecological patterns and the distribution of biomass and carbon (

The Caatinga vegetation occurs under extreme edaphoclimatic characteristics and is perhaps the most vulnerable Brazilian region to climate change (

Thus, there is a need to develop techniques to estimate and map biomass and, in turn, produce reliable information on the potential aboveground stock of forest formations and support estimates of carbon dynamics. This work examines the potential use of global climate datasets developed from satellite observations to map patterns and potential distribution of aboveground biomass in Caatinga vegetation. Unlike wetlands with high biomass driven primarily by climatic factors (

These regional characteristics favor the use of environmental variables to develop ecological and distribution models. The available remote sensing data provide spatially refined information on landscape and vegetation heterogeneity over the Caatinga SDTFs that can be readily incorporated into models to predict the potential distribution of biomass. These models are strictly mathematical or based on specific ecological theories (

The area covered by this study is located in the state of Bahia (

BSh climate according to the Köppen climate classification (

There is irregularity in the rainfall distribution in time and space. This region receives from 300 to 1180 mm of annual rainfall; however, an expressive part of the area receives average rainfall of less than 750 mm year^{-1}. Moreover, the precipitation is less than 400 mm per year in a small part. About 50% of these annual totals are generally concentrated in just three months (

Temperatures also vary in space, with annual averages between 16 and 29 °C, however, they are predominantly above 24 °C (

The soil orders that occur are: Neosols (typical in the semi-arid environment of northeastern Brazil, characterized by low natural fertility, low water retention capacity and low levels of organic matter), Latosols, Planosols, and even small fragments of Cambisols, Argisols and Luvissols (

The most representative vegetation is steppe savanna, which has a deciduous, gray and thorny physiognomy in the dry season, interspersed with cactus and bromeliads. This landscape is interrupted by seasonal semi-deciduous and deciduous forest fragments which occupy regions of soils with rocky outcrops and karst relief associated with more clayey soils, and have the presence of limestone in the valleys and slopes. The latter is more expressive on the west side of the study area. There is the occurrence of savanna and small dense ombrophilous forest fragments in mountainous regions, swamps and milder climatic pockets, common in the center-south of the state of Bahia (

The database for biomass estimation was obtained from forest surveys carried out within the study area, which include a total of 836 plots distributed in 32 forest fragments (

Data from 10 permanent plots of up to 1 ha were added to this database, registered in ForestPlots.net (

The Contendas do Sincorá National Forest was still considered as a single forest fragment in a survey carried out in 2015, inventoried by a partnership between the Caatinga Forest Management Network and the Forestry Soil and Forest Ecology and Protection Laboratories of the University State of Southwest Bahia (UESB). The measurement of arboreal individuals with DBH ≥ 15 cm was performed for this survey.

The equation developed by

(R^{2} = 92%, RMSE = 3.8 kg, bias = 7.6%) in which AGB is the aboveground dry biomass for each tree individual (kg), DBH is the diameter measured at 1.30 m above the ground (cm) and R² is the coefficient of determination. RMSE (root mean squared error) is the measure that calculates the “root mean squared error” of the errors between observed values (actual) and predictions (hypotheses).

The biomass stocks of trees per plot were summed and were extrapolated in ton per hectare (Mg ha^{-1}) according to the size of the plots. Due to the heterogeneity in the biomass stock of the plots within the same fragment, the dataset used in the analyzes was the average of the biomass stock of the plots in each of the 32 forest fragments (

The biomass stock mapping was performed using regression kriging (RK) and Inverse Squared Distance (IQD) techniques. Kriging with regression combines two processes: elaboration of a global map which illustrates the spatial behavior of the biomass, generated from a multivariate regression model; and subsequent application of ordinary kriging on the residuals generated by the regression model.

In order to introduce specific characteristics of each micro-region of the study area and provide a continuous data base, the biomass was modeled as a function of spatial variables obtained in raster layers (

After estimating the biomass stock, the data were associated with a set of predictor variables at each location. Twenty-eight geospatial covariates were selected and grouped into climatic and topographical subsets. The covariates were obtained using satellite remote sensing and globally distributed terrestrial weather stations in raster format. Each raster layer is a spatially explicit grid image, where each pixel represents the value of the described covariate.

The covariates were grouped into two distinct categories: topographical and climatic. Topographic covariates will include elevation, terrain slope, aspect (such as north and east), latitude (such as absolute latitude value), and a terrain roughness index (IRT). The climatic covariates will be composed of potential evapotranspiration, solar radiation, wind speed, cloud cover, and the set of 19 bioclimatic variables, 11 of which are derived from temperature: Average annual temperature, Average daytime interval (Monthly Average: max temp - min temp), Isothermal, Temperature Seasonality (standard deviation · 100), Maximum temperature of the hottest month, Minimum temperature of the coldest month, Annual temperature range, Average of the wettest quarter, Average of the driest quarter, Average of the warmest quarter of the year and Average of the coldest quarter of the year; and eight variables derived from Rainfall: Annual Rainfall, Rainfall in the coldest month, Rainfall in the driest month, Seasonality of Rainfall (Coefficient of variation), Rainfall in the wettest quarter, Rainfall in the driest quarter, Rainfall in the hottest quarter and Rainfall of the coldest quarter). The complete description of the topographic covariates can be consulted at http://earthenv.org/, and for the climatic covariates, more information can be obtained at https://www.worldclim.org/.

All variables acquired to fit the model, except latitude and longitude, were processed in the R software (

The biomass values used in the modeling were submitted to exploratory analysis to identify the presence of outliers and their influence on the regression assumptions. Multiple linear regression analysis was used to adjust the model, and the parameters were estimated using the least squares method. The stepwise technique was used based on the Akaike Information Criterion (AIC) with both directions in order to define which variables best explain the biomass to construct a multivariate model.

The Shapiro-Wilk and Breusch-Pagan tests were applied to verify if the model residuals presented normality and homoscedasticity, respectively. In addition, the model was submitted to the Variance Inflation Factor (VIF) test to analyze possible correlations between the explanatory variables and ensure that the model is free of multicollinearity.

The coefficient of determination (R^{2}), Mean Absolute Error (MAE) and graphical analysis of residuals were evaluated to test the accuracy of the fitted model. We also used the scatter plot of the predicted values around the 1:1 line to observe the behavior of the predictions made by the model. All analyzes were performed using the R software program (

The maps of the variables selected in the multivariate model were designed for the study area in raster format, with a spatial resolution of 1 hectare (100 × 100 m). The regression model was applied on these cells, making it possible to estimate the biomass stock in each pixel and consequently generate a global predictive map of aboveground biomass. However, this map needs to be detailed in the regression kriging procedure by adding a spatial distribution map of the regression model residuals.

Next, an exploratory analysis of the data was used to apply geostatistics on the residuals in order to know their distribution, trends and identify atypical observations (outliers). This preliminary analysis is fundamental for decision making in geostatistical procedures.

The classic estimator of

The spherical, exponential, and Gaussian theoretical models were directly fitted to the data using the Maximum Likelihood method, considering the stationarity assumption of the intrinsic hypothesis (

The theoretical model for the interpolation was selected by evaluating the Akaike Information Criterion (AIC), the degree of spatial dependence, the reduced mean error (RME) and the standard deviation of the reduced mean error (SDRME) provided by Jackknife cross-validation (

The predictive kriging quality was assessed by the mean error (ME), mean standard error (MSE), and root mean squared error (RMSE), which measure the accuracy of estimates and are provided by cross-validation by the

The prediction of biomass using the deterministic interpolator was also performed using the Inverse Weighted Distance (IWD) with exponent 2, a method known as Inverse Squared Distance (ISD). Neighborhood parameters were defined for a minimum of 10 points and a maximum of 32 points, which corresponds to the total number of sampling points in the study area. Interpolation by ISD was performed in the ArcGIS® software program (

The estimated biomass stock values in the studied fragments ranged widely, from 2.85 to 80.88 Mg ha^{-1}. There is similarity between the mean (25.17 Mg ha^{-1}) and median values (20.33 Mg ha^{-1}), which points to an approximately symmetrical data distribution, but the superiority of the mean indicates that the highest amounts of biomass are further from the center compared to the lowest stocks. The marked difference between the minimum and maximum values resulted in a high coefficient of variation (73.7%), reflecting the heterogeneity of the study area. This variation portrays how the conditions of the sampled fragments are different. Different successional stages and anthropization degrees occur in these fragments, leading to variation in the number of individuals per hectare, in the richness and diversity of species, and in the diameters and heights of the trees.

The average biomass found was similar to that estimated by ^{-1}), and lower than a fragment in Sergipe, which exhibited a contribution of 54.93 Mg ha^{-1} (^{-1} in Petrolina (PE), while ^{-1} of biomass. The discrepancy between the biomass stocks observed in the referenced works evidences the physiognomic variation that occurs along fragments of the Caatinga.

These differences in stocks between Caatinga areas are associated with precipitation, irregular rainfall distribution and the successional stage of the forest, which promote variation in the species biomass accumulation and distribution (

We identified fragments with outliers for biomass removed from the database used for subsequent analyzes as they influenced the regression assumptions. Thus, 27 fragments were used for modeling and geostatistical analysis.

The parameters used in the regression model showed significant coefficients (

The model presented a coefficient of determination (R^{2}) equal to 46% and Mean Absolute Error (MAE) of 5.9%, which indicated that the model performed well. The MAE demonstrates the model’s ability to make estimates closer to the real ones when its value tends to zero (^{2} demonstrates how much the model is able to explain the observed data. The value found herein is considered acceptable given the magnitude of the study area and the wide variation that occurs between the biomass stocks of the studied fragments. Similar conditions are also reported by ^{2} equal to 53% for a geographic model that estimates the carbon stock in the state of Minas Gerais.

The regression residuals follow a normal distribution (Shapiro-Wilk, p-value = 0.32) and are homoscedastic (Breusch-Pagan, p-value = 0.36). These characteristics are essential to define the adequacy of the model and allow adopting the Maximum Likelihood method in fitting the semivariogram, which requires

Elevation, average annual temperature and average annual precipitation are the main variables responsible for the variation in the Caatinga biomass stock in the state of Bahia (

The prediction error associated with the developed model may be due to the precision degree of applying the allometric equation used in this study, which presents R^{2} equal to 92% (

It should be noted that the accumulation of biomass in forests does not only respond to the variables considered, but also to other environmental conditions, such as soil fertility, water deficit, seasonality of rainfall and successional stages (

The semi-variance calculated in the four directions did not show significant differences and the semivariogram was considered isotropic, indicating that the spatial dependence of the regression residuals only depends on the distance between the points and is the same in all directions (

The theoretical semivariograms fitted to the data by the Maximum Likelihood method are represented in Fig. S2 (Supplementary material). It is noteworthy that the residuals had a normal distribution, which legitimizes the use of this method (

We can see the smallest nugget effect associated with the spherical model necessary to provide more accurate estimates in the kriging interpolation in Fig. S2 (

The range found was 21.000 meters and suggests that the residuals are spatially correlated up to this distance. The high range value is due to the distances between the sampled fragments and the magnitude of the study area. The relationship between the nugget effect and the threshold expresses the degree of spatial dependence of the variable, which is classified as strong (< 25%), moderate (25 to 75%) and weak (> 75%), according to

The ordinary kriging mapping of residuals by the spherical model is shown in

The error metrics associated with the interpolated surface in the ordinary kriging of residuals are presented in ^{-1} (

The combination of the regression model estimates (

There is generally an increase in the biomass stock in the north-south direction. Biomass stocks below 20 Mg ha^{-1} occur in the north of the study area due to lower precipitation rates and higher temperatures. The region surrounding the municipality of Euclides da Cunha is classified by

A small area near the municipality of Xique-Xique stood out for exhibiting higher biomass values. The junction of the residual map shows the underestimation given by the global map in this region, reinforcing that the combination between a global interpolator and a geostatistician better represent the spatial distribution of a variable. The evident detail given in this specific area is due to the higher concentration of fragments sampled there, which in turn brought the regression kriging map (

The interpolated surface presented mean error (ME) values of -0.38, which suggests an overestimation of the biomass stock prediction by the ISD. The ME measures the model’s tendency to under or overestimate a variable of interest and should be close to zero for unbiased predictions (^{-1}. Although this type of interpolation is useful in the absence of spatial structure of the variable of interest, or even in the presence of weak dependence structure, the use of the ISD was relevant because it showed that the inconsistency in the mapping by kriging with regression was influenced by the subsampling in the study area.

Small isolated areas with smaller stocks can be observed in the south-central portion of the map where biomass stocks are greater than 25 Mg ha^{-1}, which constitutes an effect of the low estimated biomass in the fragments sampled there. This shows that these fragments do not capture the real variation of biomass stock that occurs in the study area and reinforces that the underestimation of the Chapada Diamantina represented in the regression kriging map is due to the effect of these fragments. Despite this, in visually comparing the maps it is notable that the regression kriging is able to capture more specific details regarding the biomass distribution in the studied area than the ISD, which makes it more efficient in the mapping.

The largest biomass stocks are in the south of the map. The biomass in the municipality of Contendas do Sincorá was estimated at 37 Mg ha^{-1}, while ^{-1}. The authors classify the predominant vegetation in this area as arboreal Caatinga and in a late successional stage. Considering that these authors used the direct quantification method which is more accurate, it is stated that the biomass in the present study was slightly overestimated in this region.

The highest biomass stocks indicated in southern Bahia suggest that these are potential areas in the Caatinga to store carbon. Considering that the carbon contents in the biomass are generally around 50% (

Low biomass values were found for the Chapada Diamantina region, where the municipalities of Morro do Chapéu and Lençóis are located. Chapada Diamantina is characterized by a set of mountains that reach more than 1000 meters in elevation, with high precipitation levels and milder temperatures (^{-1} in this region and classified it as one of the areas with the highest plant biomass in the entire Caatinga.

However, this region does not support a homogeneity of dense vegetation. Chapada Diamantina is a set of communities which form a mosaic rich in physiognomies, depending on the topography, nature and depth of the soil. The high altitudes and busy relief result in a high erosion and occupation rate of the areas by rocky outcrops, which is limiting to vegetation development, with formations having different physiognomies from the general dominant context of this region occurring in these soils (

Therefore, it is likely that the low biomass values in this area on the map (

However, lower mapping accuracy was expected due to the richness of physiognomies of the Caatinga in Bahia, as characterizing the biomass in heterogeneous environments such as tropical forests is a challenge (

There is no consolidated database available which covers the entire territory of the state of Bahia, which validates the data compiled in this study. Although not as consistent, the data provide the potential for biomass in each patch. A base of inventories capable of capturing the entire spectrum of vegetation variation in Bahia and showing its plurality of biomes, successional stages, anthropization degrees and phytophysiognomies is necessary. The existence of a consistent basis would enable developing techniques and studies capable of supporting the conservation and sustainable exploitation of the state’s forest resources with greater precision.

This information would provide greater detail in the final map in the present study and would promote a substantial improvement in mapping.

Due to the heterogeneity of the Caatinga vegetation in Bahia, it is suggested to consider the different phytophysiognomies, ecoregions or climatic zones of the biome as independent variables in the biomass modeling in future works. This distinction has been used to reduce the spatial variation common in vegetation data over large areas (

Similarly,

Our analysis explored a variety of data sources and analytical tools that can be applied to develop predictive maps that incorporate the observed spatial pattern of biomass in landscapes at the scale of the Caatinga biome in the state of Bahia. However, the context of available data still lacks more excellent coverage or prescription for some areas, therefore some limitations remain.

Firstly, it is common to face problems using inventory data collected differently in each location, mainly associated with taxonomically reliable data and correctly defined technical protocols in a given location. As a result, some trees may be measured differently, and individual biomass estimates may show under- or over-estimate trends. For example, many local, regional, and pan-tropical equations use the diameter at 1.30 m ground level (dbh) in biomass estimates, and many inventories only have a diameter at 0.30 m ground level (diameter at base). Furthermore, not all Caatinga forest fragments or land use forms can be equally represented, leading to high standard deviations or region-biased biomass estimates. Although the locations were spatially well distributed for the Caatinga biome in Bahia, and some locations were close to each other, an effort was made to correct and standardize dendrometric measurements in specialized protocols and to use forest land-use maps at a scale of finer detail. The analysis of the biomass stocks discussed is also based on accurate classifications of land cover and markedly anthropized gradients, which may not be available or may be misaligned with the inventoried plots. This helps to correct errors about the representativeness and potential distribution of biomass at the plot level and the biome scale.

Second, different scales of environmental data (10 × 10 km tables) and inventory form another complication, as plots covering a few hectares are unlikely to represent an area of 10 × 10 km (

Thirdly, and related to the available environmental layers and their potential to explain the biomass in Caatinga vegetation, it was noted that the functional importance of many environmental variables suggests poor results among all the variables tested based on characteristics related to predictive importance. The high collinearity may help to explain this fact; however, other modeling approaches, as well as other environmental variables, should be tested, mainly related to the soil of the biome. For example, edaphic factors are vital in explaining many attributes and ecological patterns of the forest, mainly related to productivity and diversity in tropical forests (

We integrate

The quality of the regression model suggests that three main environmental factors possibly govern aboveground biomass. Increases in temperature and altitude suggest a reduction in biomass. However, the combination of average annual precipitation and moderate altitude conditions is the ideal climate to support the highest biomass stocks. Additional studies with a larger sample population and other variables can improve the model. Compared to ISD, regression kriging revealed in more detail the spatial variability of biomass. The adoption of a higher sampling intensity has the potential to maximize this detail given by the technique. A robust base of forest inventories is needed for Bahia in order to enable studies that allow for a more accurate understanding of the potential of the state’s forest resources.

Finally, we mapped the potential biomass distribution of the Caatinga biome in Bahia, but we do not guarantee that the conservation or preservation of these is established. The overwhelming effect of human activities (deforestation, fires, and the advance of agriculture and livestock) has led to substantial reductions in biomass and carbon stocks. Given the enormous carbon-storage capacity of tropical dry forests, the next avenue for research is to model the effects of climate change and future land-use changes on biomass patterns to refine conservation prospects in a changing world. These efforts are urgently needed in light of the constant and increasing rates of deforestation in the Caatinga biome in recent years and potentially threaten the planet’s high biodiversity hotspots.

We thank the ForestPlots.net, Inema and Caatinga forest management network teams for permission and availability of data to conduct the research. This project was supported by ForestPlots.net approved Research Project #142: “Mapping the biomass and carbon stock, volume and diversity in tropical forests", and the UK NERC Newton Fund project “Nordeste” (NE/N012550/1). The authors thank the Coordination for the Improvement of Higher Education Personnel (Capes) and the National Council for Scientific and Technological Development (CNPq) for financial support in a research grant for the first and third authors of this work.

The conceptualization, data curation and formal analysis: Santos HKV, Lima RB. The Investigation, Methodology, Validation, Visualization and Project administration: Santos HKV, Lima RB, Silva TTS, Souza RLT, Oliveira CP, de Paula A, Barreto-Garcia PAB, Veenendaal E, Queiroz L, Moonlight P, Rodrigues PMS, Santos RB, Sarkinen T, Pennington T, Domingues T, Cardoso D, Phillips O. Writing-original draft: Viana HK, Lima RB. Writing review and editing: Viana HK, Lima RB, Moonlight P, Phillips O, Cardoso D.

Location of the sampled fragments in the Caatinga biome, in the state of Bahia, Brazil.

Flowchart of Regression Kriging (GWR-K) and AGB mapping predictions in this study.

Mapping by ordinary kriging of the residuals of the regression model.

Variables selected by the multivariate regression model to estimate biomass in the Caatinga, Bahia.

Global map of biomass stock obtained by the regression model (A), interpolated by kriging with regression - RK (B) and interpolated by the Inverse Square of Distance technique - IQD (C).

Bioclimatic variables and elevation above sea used to adjust the biomass prediction model in the caatinga, Bahia.

Variables | Descriptions |
---|---|

BIO 1 | Average annual temperature (°C) |

BIO 2 | Average daytime range (monthly average) (°C) |

BIO 3 | Isothermal (Bio 2/ Bio 7) (·100) (°C) |

BIO 4 | Temperature seasonality (standard deviation·100) (°C) |

BIO 5 | Maximum temperature of the warmest month (°C) |

BIO 6 | Minimum temperature of the coldest month (°C) |

BIO 7 | Override temperature range (Bio 5 - Bio 6) (°C) |

BIO 8 | Average temperature of the wettest quarter (°C) |

BIO 9 | Average temperature of the driest quarter (°C) |

BIO 10 | Average temperature of the warmest quarter (°C) |

BIO 11 | Average temperature of the coldest quarter (°C) |

BIO 12 | Average annual rainfall (mm) |

BIO 13 | Rainfall in the wettest month (mm) |

BIO 14 | Rainfall of the driest month (mm) |

BIO 15 | Seasonality of precipitation (Coefficient of variation) (mm) |

BIO 16 | Precipitation of the wettest quarter (mm) |

BIO 17 | Rainfall in the driest quarter (mm) |

BIO 18 | Precipitation of the warmest quarter (mm) |

BIO 19 | Rainfall of the coldest quarter (mm) |

Elev | Elevation (m) |

Coefficients estimated by the regression model and Variance Inflation Factor (VIF).

Parameter | Variable | Coefficient | p-value (<0.05) | VIF |
---|---|---|---|---|

b_{0} |
intercept | 519.3 | 0.000346 | - |

b_{1} |
Elev | -0.0000594 | 0.000293 | 1.70 |

b_{2} |
Temp | 1.831 | 0.045831 | 1.18 |

b_{3} |
Prec | -0.03326 | 0.016238 | 1.59 |

Statistical parameters and criteria of the theoretical models adjusted for the residuals of the regression model. (C_{0}): nugget effect; (C_{1}): threshold; (A): range; (GD): degree of spatial dependence; (EMR): mean reduced error (S_{ER}): standard deviation of reduced errors (AIC): Akaike Information Criterion.

Model | C_{0} |
C_{1} |
A (m) | GD (%) | EMR | S_{ER} |
AIC |
---|---|---|---|---|---|---|---|

Exponential | 11.75 | 49.00 | 21000 | 23.90 | 0.018 | 1.001 | 187.3 |

Gaussian | 13.69 | 48.12 | 21000 | 28.45 | 0.020 | 1.002 | 186.3 |

Spherical | 5.06 | 45.45 | 21000 | 11.14 | 0.016 | 0.992 | 185.4 |

Statistics of the interpolation by ordinary kriging of the residuals of the regression model. (ME): Mean Error; (MSE): Mean Standardized Error; (RMSSE): Root Mean Square Error Standardized; (RMSE): Root Mean Square Error; (ASE): Mean Standard Error. For the semivariogram model to be suitable for interpolation, ME and MSE must be close to 0; the RMSSE should be close to 1 e; the RMSE and the ASE should have similar values and the smallest possible (

Statistics | Value |
---|---|

ME | -0.16 |

MSE | 0.01 |

RMSSE | 0.97 |

RMSE | 6.39 |

ASE | 6.82 |

Fig. S1 - Graphical distribution of regression residuals (a); and observed

Fig. S2 - Fitting of the theoretical semivariogram models for the residuals of the regression model, using the Maximum Likelihood method.