Landsat-8 and Sentinel-2 Based Prediction of Forest Plantation C Stock Using Spatially Varying Coefficient Bayesian Hierarchical Models

Abstract

This study sought to establish the performance of Spatially Varying Coefficient (SVC) Bayesian Hierarchical models using Landsat-8, and Sentinel-2 derived auxiliary information in predicting plantation forest carbon (C) stock in the eastern highlands of Zimbabwe. The development and implementation of Zimbabwe’s land reform program undertaken in the year 2000 and the subsequent redistribution and resizing of large-scale land holdings are hypothesized to have created heterogeneity in aboveground forest biomass in plantation ecosystems. The Bayesian hierarchical framework, accommodating residual spatial dependence and non-stationarity of model predictors, was evaluated. Firstly, SVC models utilizing Normalized Difference Vegetation Index (NDVI), Soil Adjusted Vegetation Index (SAVI), and Enhanced Vegetation Index (EVI), derived from Landsat-8 and Sentinel-2 data and 191 sampled C stock observations were constructed. The SVC models built for each of the two multispectral remote sensing data sets were assessed based on the goodness of fit criterion as well as the predictive performance using a 10-fold cross-validation technique. The introduction of spatial random effects in the form of Landsat-8 and Sentinel-2 derived covariates to the model intercept improved the model fit and predictive performance where residual spatial dependence was dominant. For the Landsat-8 C stock predictive model, the RMSPE for the non-spatial, Spatially Varying Intercept (SVI) and SVC models were 8 MgCha⁻¹, 7.77 MgCha⁻¹, and 6.42 MgCha⁻¹ whilst it was 7.85 MgCha⁻¹, 7.69 MgCha⁻¹ and 6.23 MgCha⁻¹ for the Sentinel-2 C stock predictive models, respectively. Overall, the Sentinel-2-based SVC model was preferred for predicting C stock in plantation forest ecosystems as its model provided marginally tighter credible intervals, [1.17–1.60] MgCha⁻¹ when compared to the Landsat-8 based SVC model with 95% credible intervals of [1.13–1.62] Mg Cha⁻¹. The built SVC models provided an understanding of the performance of the multispectral remote sensing derived predictors for modeling C stock and thus provided an essential foundation required for updating the current carbon forest plantation databases.

1. Introduction

Since the onset of the Fast Track Land Reform (FTLRP) program in Zimbabwe and the subsequent redistribution and resizing of large-scale land holdings in the year 2000, plantation forests within and around the neighborhood of resettlement areas continuously faced physical distress [1]. Zimbabwe lost close to 224,000 ha of tree cover, which is equivalent to a forest degradation of 17% and 88 Mt of CO₂ emissions from 2000 to 2021 [2]. The steady decline of land under forest over the years has therefore been mainly driven by the activities of the year 2000 taking place on resettled land, which included wildfires, illegal logging, mining, and agriculture expansion. Yet the amount of additional biomass that can be accumulated in these areas depends much on the forest condition and land management practices.

A number of studies modeling the relationship between C stock and remote sensing-derived information adopt a standard approach where the effects of independent variables, including vegetation indices on C stock, remain constant over space in the model [3,4]. Given the history of Zimbabwe’s land reform program, it is well known that management and conservation practices of plantation forests have been changing over the years, particularly since the adoption of the program in the year 2000. Greater portions of formally designated commercial plantation forests were occupied and replaced by subsistence farmers across the country. It is upon this premise that it cannot be assumed that the impact of covariates on C stock is constant in these ecosystems. The carbon sequestration potential in various regions of the occupied forest plantations, therefore, remains unknown as the management and conservation practices of the managed ecosystems have been largely modified and altered [5]. Yet government and other stakeholders in the timber industry are obliged under the 2015 Paris climate agreement to provide accurate accounts of the country’s carbon sequestration potential for managed and natural forest ecosystems.

The literature proffers methods used in the estimation of AGB, either as direct or indirect approaches from forest inventory data. These methods either use allometric equations, conversion factors such as wood density, or biomass expansion factors (BEF) [1]. In spite of the advantages of using conventional methods given in [6,7] as generally providing accurate AGB estimates, these techniques are also regarded as time-consuming and environmentally unfriendly. The inaccessibility of some areas due to complicated topography and forest conditions also makes conventional methods less attractive to AGB estimation, especially in extensive areas [8]. However, remote sensing is emerging as a promising technique free of the aforementioned limitations as it offers cost-effective methods of AGB monitoring through stratification of canopy density and forest types. Repeated application of remote sensing leads to the generation of historical data critical for change detection analysis and incorporation into a Geographic Information System (GIS) for integration with other datasets [9]. However, remote sensing methods are also not immune from limitations as they usually face limitations in areas of bad weather conditions, especially in areas with cloud cover, and require a certain level of training for effective use and application [7].

Remote sensing C stock data derived from satellite imagery has significantly grown over the years. This data has been the basis for informing international policy agreements associated with CO₂ emissions into the atmosphere, mainly from deforestation and other land use land cover changes (LULCC) [10]. Remote sensing is regarded as a powerful tool for deriving AGB and forest structure as it offers practical means of acquiring spatially-distributed forest carbon from local, continental, and global scales [5,11]. Forest biomass can generally be measured from three broad categories of remote sensing (RS) data, namely, passive optical remote sensing, radio detection and ranging, microwave (radar), and light detection and ranging (LiDAR) [8,10]. Passive optical spectral reflectance is responsive to vegetation structural attributes (tree density, leaf area index, and crown size), shadow, and texture [8,12]. Crown size, tree density, and leaf area index (LAI) are strongly correlated with AGB. Radar remote sensing measure geometrical and dielectric attributes of forests. LiDAR RS methods of biomass measurement characterise forest vertical structure and height. Remote sensing (RS) and ancillary technologies such as Geographical Information System (GIS) are practical and cost-time effective, allowing for imaging and research on extensive and inaccessible areas.

The estimation of plantation forest C stock, together with other structural parameters using new regeneration multispectral remote sensing data, is a relatively new ground in the climate change and carbon monitoring arena. Research in this field has demonstrated the potential for remote sensing data as tools for developing estimates of forest attributes, amongst them, being Above Ground Biomass (AGB), either as a standalone or coupled with other earth observation techniques [11,13]. Spatial regression models applied in mapping C or biomass using new-generation multispectral remote sensing data may fail to clearly make room for residual spatial dependence [14,15]. Modeling of natural resource variables without explicitly accommodating spatial variation can be justified if covariates can account for all the spatially structured dependence. Such assumptions do not hold in many practical applications involving georeferenced data.

Failure to account for spatial dependence in the modeling of regionalized variables can lead to inaccurate model parameters and incorrect predictions [9]. Subsequently, disregarding spatial correlation in electing model choices can result in higher prediction uncertainty for inference of the outcome variable. In addition to the drawbacks highlighted in the foregoing, non-Bayesian spatial modeling can further lead to underestimation of uncertainty [9,16] as traditional spatial regression estimation methods assume stationarity of the covariance matrix, $\mathsf{\Sigma }$. The ultimate effect is the derivation of standard error estimates that are unable to take all the uncertainties in the parameters into account. Checking for spatially correlated residuals when spatial data are employed in the modeling of above-ground biomass is therefore critical.

Some attention has been given to spatially varying coefficient (SVC) models in the literature [17,18,19]. The Bayesian framework of statistical inference is the bedrock of SVC models in which analyses make use of samples derived from Markov Chain Monte Carlo (MCMC) techniques from the posterior distribution of model parameters [20]. What makes Spatially Varying Intercept (SVI) and SVC models unique in applied geostatistical and the remote sensing literature is their ability to consider the residual spatial dependence and the non-homogeneity in model parameters differently than ordinary geostatistical approaches. SVI models assume the model intercept is spatially varying, whilst the SVC models assume all the model regression coefficients to be spatially varying [21]. A number of applications and methodologies utilizing spatially varying coefficient models are documented in the literature. Amongst them are the geographically-weighted regression (GWR) by [21], who employed the technique of canopy height prediction using remote sensing data. Application of these models resulted in significant improvement in model fit. On the other hand, [22] made provisions for spatial dependence and parameter non-homogeneity by exploring geostatistical kriging variants for forest canopy height prediction. Co-kriging and regression kriging models resulted in significant improvements in model fit.

However, it has been shown in recent times that GWR might not be robust to correlation among predictors and can potentially lead to inaccurate results when complex correlation structures are involved [18,23]. Again, from an inferential viewpoint, GWR can present problems when drawing inferences regarding prediction uncertainty and model parameters. The lack of valid underlying probability models in GWR makes prediction uncertainty, and standard parameter error estimates difficult to justify. For instance, uncertainty maps produced from kriging and co-kriging techniques make no consideration of the uncertainty in the variogram-based spatial covariance parameters. This is an established and common weakness encountered when using these geostatistical approaches [24]. It is possible to estimate spatial covariance parameters within the SVC and SVI models using a Bayesian hierarchical construction. Such an approach allows the propagation of uncertainty to the prediction of the response variable [25]. In such scenarios, a better and statistically defendable map of uncertainty can be produced than would else be produced when GWR or traditional kriging methods are utilized.

Landsat-8 and Sentinel-2 are multispectral platforms categorized as new generation remote sensing sensors with enhanced spectral and spatial properties than previous missions of the Landsat series. It is this perceived improvement in earth observation properties that we expect to give a dividend to C stock predictive models for applications in carbon accounting and inventorying under the United Nations Framework Convention on Climate Change (UNFCC) [8]. We employed a Bayesian hierarchical framework with spatially varying coefficients (SVC) using predictor information derived from the aforementioned sensors for predicting C stock. The modeling framework considered the non-stationarity and residual spatial dependence of model covariates through the inclusion of spatial random effects into the SVC models. We, therefore, developed and compared the performance of C stock predictive models under spatially varying regression coefficients derived from Landsat-8 and Sentinel-2 predictors. Prediction accuracy and uncertainty quantification for the REDD collaborative program in the developing world is a critical aspect of the Carbon Measurement, Reporting, and Verification Systems (MRVs) of the United Nations (UN-REDD, 2009; CMS, 2014). Thus, in this paper, we explored how spatially varying coefficient models constructed using a Bayesian hierarchical set-up with aiding information from Landsat-8 and Sentinel-2 multispectral sensors and implemented through Markov Chain Monte Carlo (MCMC), perform in C stock prediction in plantation forest ecosystems.

2. Materials and Methods

2.1. Study Area

The study was carried out at lot 75 A of Nyanga Downs in Nyanga district in the Eastern Highlands of Zimbabwe (Figure 1). The study area is dominated by Eucalyptus grandis, Eucalyptus camaldulensis, and Pinus patula plantation forest species which have some of its patches being used for agriculture, grazing, and gold panning and is located between latitude 32°40′E and 32°54′E and 18°10′12″S and 18°25′4″S longitude as illustrated in Figure 1 [5,26]. Grazing, agriculture, and gold panning activities came after part of the commercially owned plantation forests were redistributed to small and medium-sized indigenous farmers in 2000. This development has increased the interface between settlements and timber plantations in all forests originally designated under forest plantations in Zimbabwe. The study area covers an area of approximately 2766 ha. Rainfall amounts are varied, with a mean annual precipitation range of 741 mm to 2997 mm. Annual mean temperatures range from a minimum of 9 °C to 12 °C to a maximum of 25 °C to 28 °C. The weather is very hot, and extensive wildfires occur in the high-elevation grasslands from August to November when the grasses are dry [26].

As illustrated in Figure 1c, Eucalyptus camaldulensis and Pinus patula are the most dominant species in the study area. Pockets of cultivated land within the plantations are evident and are partly responsible for the present biomass density in the sampled region. Greater portions of former plantation vegetation have been cleared by pockets of resettled small-scale and medium-scale farmers venturing into coffee and tea plantations and, in some cases, for dairy farming. Patches of unprotected forest plantations are still present but remain vulnerable to attack for agriculture by settled farmers in the area.

2.2. Sampling Design

We employed the spatial coverage sampling scheme that exploits the Mean Squared Shortest Distance (MSSD) for the optimization of data locations in the study domain. The k-means clustering scheme for equal area coverage was therefore used [27] for obtaining a representative sample from the studied region. The work of [28] demonstrated how the mapping and estimation of mean spatial problems could be resolved through the employment of a uniform coverage sampling scheme. MSSD is particularly suited for areas where sampling campaigns cannot be extended beyond a single phase. As illustrated in Figure 2, the smallest separation distance between samples was approximately 8 m whilst the largest separation distance between samples was 2500 m.

2.3. Above-Ground Biomass Measurements and C Stock Derivation

We sampled and collected measurements of Diameter at Breast Height (DBH) for Above Ground Biomass (AGB) estimation for all trees with a DBH of more than 10 cm (1.3 m) using 500 m² circular plots from the 19 September 2021 to 24 October 2021 in Manicaland Province of Zimbabwe as illustrated in Figure 1. The sampling program meant to measure DBH for subsequent AGB and C stock estimation utilized the MSSD optimization function, resulting in 191 sampling plots being measured from the study area. The 191 sample plot measurements of DBH were then transformed into per plot C stock data using allometric equations of [4] for the Pinus species, whilst the allometric equations of [10] were used for deriving C stock for the Eucalyptus species. The aforementioned allometric equations were also applied by [26].

Allometric equations shown in Equations (1) and (2) were used for the calculation of Above Ground Biomass (AGB) for the Pinus patula and the Eucalyptus grandis and Eucalyptus camaldulensis species, respectively [1]. A default conversion factor of 0.47 used by the IPCC was applied to derive AGB to C stock.

$$tDw=e^{\left(-1.170+2.119\times ln\left(dbh\right)\right)}$$

(1)

$$tDw=0.39\times {\left(dbh\right)}^{2.142}$$

(2)

2.4. Modelling Framework

It is a common geostatistical practice to assume at location $S\in D\subseteq {ℝ}^2$ where $s$ is a vector of observed $x,y$ coordinates within the domain $D$. A Gaussian response variable $y\left(s\right)$ is therefore modeled through the regression model as in Equation (3):

$$y\left(\mathit{s}\right)=x{(\mathit{s})}^{\prime }\beta +\tilde{\mathit{x}}{(\mathit{s})}^{\prime }\mathit{w}\left(\mathit{s}\right)+\epsilon \left(\mathit{s}\right)$$

(3)

$x{(\mathit{s})}^{\prime }$ denotes a set of $p$ covariates in the model. In this case, the linear mean structure accounting for wide-scale variation in the response is comprised of $px1$ vectors of $\mathit{x}\left(\mathit{s}\right)$ which include an intercept and spatially varying georeferenced predictors together with an associated column vector of model coefficients $\mathit{\beta }={({\beta }_0,{\beta }_1,\dots ,{\beta }_{p-1})}^{\prime }$. The $\tilde{\mathit{x}}\left(\mathit{s}\right)$ in the model represents a $q x 1$ vector accommodating the intercept and those predictors from $x\left(s\right)$ whose impact on the response is assumed to vary spatially. The space-varying impact is obtained from the vector of spatial random effects $\mathit{w}\left(\mathit{s}\right)={\left(w_1\left(\mathit{s}\right),w_2\left(\mathit{s}\right),\dots ,w_q\left(\mathit{s}\right)\right)}^{\prime }$. The specification of different combinations of $\tilde{\mathit{x}}\left(\mathit{s}\right)$ and the associated $\mathit{w}\left(\mathit{s}\right)$ leads to the formation of different sub-models. We model $\epsilon \left(\mathit{s}\right)$ as an independent white noise process that takes care of the micro-scale (measurement error) variation. As such, with the collection of $n$ C stock locations, $S={\mathit{s}}_{\mathit{1}},{\mathit{s}}_{\mathit{2}},\dots ,{\mathit{s}}_{\mathit{n}}$, we assume the $\epsilon \left({\mathit{s}}_{\mathit{i}}\right)$’s are independent and identically distributed ($iid$) as provided by $N\left(0,{\tau }^2\right)$ where ${\tau }^2$ is the nugget.

The spatial structure of this model is generally introduced by way of a multivariate Gaussian process (GP) [24,25,29] in which a cross-covariance function clearly models the covariance of $\mathit{w}\left(\mathit{s}\right)$ within and among data points. The added flexibility in this model is documented in the literature [17,30]. We assume in this study that elements of $\mathit{w}\left(\mathit{s}\right)$ emanate from $q$ independent univariate GPs. Precisely, the process associated with the $k-\mathrm{th}$ prediction is $w_k\left(\mathit{s}\right)~GP\left(0,\left(.,.;{\mathit{\theta }}_k\right)\right)$ where $C\left(\mathit{s},{\mathit{s}}^{*};{\mathit{\theta }}_k\right)=Cov\left(w_k\left(\mathit{s}\right),w_k\left({\mathit{s}}^{*}\right)\right)$ is a valid spatial covariance function modelling the covariance related to a pair of observations $\mathit{s}$ and ${\mathit{s}}^{*}$. The process outcomes are gathered into an $n x 1$ vector, say ${\mathit{w}}_k={\left(w_k\left({\mathit{s}}_{\mathbf{1}}\right),\dots ,w_k\left({\mathit{s}}_n\right)\right)}^{\prime }$ which permits a multivariate normal distribution $MVN\left(0,{\Sigma }_k\right)$, in which ${\Sigma }_k$ is an $n x n$ covariance matrix of ${\mathit{w}}_k$ with the $\left(i,j\right)-\mathrm{th}$ element provided by $C\left({\mathit{s}}_i,{\mathit{s}}_j;{\mathit{\theta }}_k\right)$. Evidently, $C\left(\mathit{s},{\mathit{s}}^{*};{\mathit{\theta }}_k\right)$ cannot just be a function, but guarantees that the resulting ${\Sigma }_k$ matrix is positive definite and symmetric. Functions of this type are regarded as positive definite and are referred to as the characteristic function of a symmetric stochastic variable [24,25,31].

We denote $C\left(\mathit{s},{\mathit{s}}^{*};{\mathit{\theta }}_k\right)={\sigma }_k^2\rho \left(\mathit{s},{\mathit{s}}^{*};{\theta }_k\right)$ with ${\mathit{\theta }}_k=\left\{{\sigma }_k^2,{\varphi }_k\right\},\rho \left(.;{\varphi }_k\right)$ to be a valid spatial correlation function in which ${\varphi }_k$ computes the correlation decay rate and $var\left(w_k\right)={\sigma }_k^2$. We assumed for all the accompanying analyses an exponential correlation function $\rho \left(\left|\left|\mathit{s}-{\mathit{s}}^{*}\right|\right|;{\varphi }_k\right)=exp\left(-{\varphi }_k\left|\left|\mathit{s}-{\mathit{s}}^{*}\right|\right|\right)$, where $\left|\left|\mathit{s}-{\mathit{s}}^{*}\right|\right|$ represents the Euclidean distance between location $\mathit{s}$ and location ${\mathit{s}}^{*}$. Completion of the Bayesian modeling framework requires specification and assignment of prior distributions to the parameters of the model, where inference proceeds by sampling from the posterior distribution of the modeled parameters. As standard practice, we assume $\beta ~MVN\left({\mathit{\mu }}_{\beta },{\Sigma }_{\beta }\right)$ prior where ${\mathit{\mu }}_{\mathit{\beta }}=0$ and ${\Sigma }_{\beta }=\mathrm{10,000}{\mathit{I}}_p$, whilst the spatial variance components ${\sigma }_k^2$’s and the measurement error variance ${\tau }^2$ are designated inverse-Gamma, $IG\left(a,b\right)$ priors. The spatial decay parameters ${\varphi }_k~Unif\left(a,b\right)$ with the lower and upper bounds of the distribution covering the geographic domain of the sampled study region.

Applying notation similar to the ones used by [32], we can specify the model parameter posterior distribution as $p(\Omega |\mathit{y})$ where:

$$\Omega =\left\{\mathit{\beta },{\mathit{w}}_{\mathbf{1}},{\mathit{w}}_{\mathbf{2}},\dots ,{\mathit{w}}_{\mathit{q}},{\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_q^2,{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_q,{\tau }^2\right\}$$

$\mathit{y}={\left(y\left(s_1\right),\dots y\left(s_n\right)\right)}^{\prime }$ and is proportional to:

$${\prod }_{k=1}^{q}Unif\left({\varphi }_k|a_{\varphi k},b_{\varphi k}\right)x{{\displaystyle \prod }}_{k=1}^{q}IG\left({\sigma }_k^2|a_{\sigma k},b_{\sigma k}\right) x N\left(\mathit{\beta }|{\mathit{\mu }}_{\beta }\right)xIG\left({\tau }^2|a_{\tau },b_{\tau }\right) x$$

$${{\displaystyle \prod }}_{k=1}^{q}N\left({\mathit{w}}_k|\mathbf{0},{\Sigma }_k\right)x{{\displaystyle \prod }}_{i=1}^{n}N(y\left({\mathit{s}}_i\right)|\mathit{x}{({\mathit{s}}_i)}^{\prime }\mathit{\beta }+\tilde{\mathit{x}}\left(\mathit{s}{)}^{\prime }\mathit{w}\left(\mathit{s}\right),{\tau }^2\right)$$

(4)

An effective Markov Chain Monte Carlo (MCMC) function for the estimation of Equation (4) is derived through updating of $\mathit{\beta }$ from its full conditional and utilizing Metropolis procedures for the remainder of the parameters. Reparameterization of the model is an alternative way of ensuring the spatial random effects $\mathit{w}$ do not require direct sampling [33]. The spatial random effects could represent other independent variables that are spatially structured and not taken into consideration in the current modeling approach. Nonetheless, the MCMC process produces posterior samples of the parameter space, $\Omega$.

From a prediction point of view, if $S_0=\left\{{\mathit{s}}_{0,1},{\mathit{s}}_{0,2},\dots ,{\mathit{s}}_{0,m}\right\}$ is a set of $m$ new sites, the spatial random effects posterior predictive distribution corresponding to the $k-th$ regression coefficient is provided by:

$$p({\mathit{w}}_{k,0}|\mathit{y})\propto \int p({\mathit{w}}_{k,0}|{\mathit{w}}_k,\Omega ,\mathit{y})p({\mathit{w}}_k|\Omega ,\mathit{y})p(\Omega |\mathit{y})d\Omega {\mathit{w}}_k$$

(5)

where ${\mathit{w}}_{k,0}={\left(w_k\left({\mathit{s}}_{0,1}\right),w_k\left({\mathit{s}}_{0,2}\right),\dots ,w_k\left({\mathit{s}}_{0,m}\right)\right)}^{\prime }$.

Since we are making use of MCMC sample-based inference, the integral in Equation (5) does not need to be evaluated precisely. Instead, given $L$ posterior samples for the parameter space ($\Omega$), that is, ${\left\{{\Omega }^{\left(l\right)}\right\}}_{l=1}^L$, composition sampling can be used to derive this distribution [33] by first drawing $w_k^l$ followed by ${\mathit{w}}_{k,0}^{\left(l\right)}$ for each $l$ from $p({\mathit{w}}_{k,0}|{\mathit{w}}_k^{\left(l\right)},{\Omega }^{\left(l\right)},\mathit{y})$. The last distribution is multivariate normal as it is a derivative of a conditional distribution from a multivariate normal distribution. More specifically, the process realizations over the new sites are conditionally independent of the measured outcomes given the values over the sampled locations and the process parameters. Expressed differently, $p(w_{k,0}|w_k,\Omega ,y)=p(w_{k,0}|w_k,\Omega )$ is a multivariate distribution with mean and variance given by

$$E[{\mathit{w}}_{k,0}|{\mathit{w}}_k,\Omega ]=Cov\left({\mathit{w}}_{k,0},{\mathit{w}}_k\right)va{r}^{-1}\left({\mathit{w}}_k\right){\mathit{w}}_k=R_0{({\varphi }_k)}^{\prime }R{({\varphi }_k)}^{-1}{\mathit{w}}_k$$

(6)

and

$$var[{\mathit{w}}_{k,0}|{\mathit{w}}_k,\Omega ]={\sigma }_k^2\{R\left({\varphi }_k\right)-R_0{({\varphi }_k)}^{\prime }R{({\varphi }_k)}^{-1}{R}_0\left({\varphi }_k\right)\}$$

(7)

where $R_0({\varphi }_k$ is an $n x n$ matrix with $\left(i,j\right)-\mathrm{th}$ element specified as $\rho \left({\mathit{s}}_{0,i},{\mathit{s}}_j;{\theta }_k\right)$ and $R({\varphi }_k$ is an n x n matrix with $\left(i,j\right)-\mathrm{th}$ element provided by $\rho \left({\mathit{s}}_i,{\mathit{s}}_j;{\varphi }_k\right)$. Repetition of the same procedure results in the generation of samples for all the ${\mathit{w}}_{k,0}$’s. Lastly, for a set of predictors at unsampled locations ${\mathit{s}}_0$, posterior predictive distribution samples regarding the outcome variable $y{({\mathit{s}}_0)}^l$, are derived from $N(\mathit{x}{(\mathit{s})}_0^{\prime }{\mathsf{\beta }}^{\left(l\right)}+\tilde{\mathit{x}}({\mathit{s}}_0)\prime {\mathit{w}}_0^{\left(l\right)},{\tau }^{2\left(l\right)})$ for $l=1,2,\dots ,L$.

We evaluated 95% credible interval widths ($CIWs$) using the posterior predictive distribution of Landsat-8 and Sentinel-2 C stock models by calculating the difference between the 2.5% and the 97.5% quantile bounds. The 95% $CIWs$ were therefore used as summaries of the C stock prediction uncertainty for the Landsat-8 and Sentinel-2 C stock-based spatially varying coefficient models.

2.5. Competing Models

We derived five candidate models for each of the two multispectral remote sensing-based C stock models using Equation (1) using $NDVI$ and $SAVI$ as predictors of C stock [34]. The models included the non-spatial where $w_k$’s is fixed to zero; the spatially varying intercept ($SVI$) in which we only included the spatial random effects related to the model intercept; the complete $SVC$ model in which all predictors have associated spatial random effects; the $SVC-{\beta }_1$ which included the spatial random effects for the intercept and $NDVI$ predictor variables and the $SVC-{\beta }_2$ which included the spatial random effects for the intercept.

We utilized empirical semivariograms modeled for the residuals of the independent error model as guidelines for candidate model $IG$ and $Unif$ hyperprior specifications. Precisely, for the variance parameters, the Inverse Gamma hyperprior $a$ was set equal to $1.76$, which would result in a mean prior distribution equal to $b$ and infinite variance [35]. To add on, the models’ $b$ hyperpriors for the ${\tau }^2$ and ${\sigma }^2$’s were calibrated in accordance with the sample variograms of the simple linear regression model residuals derived from Landsat-8 and Sentinel-2 sensor’s nugget and partial sill. We programmed the prior for the spatial decay parameter $\varphi$’s to $Unif \left(0.38, 0.0012\right)$ which, adopting the exponential covariance function, equates to support an effective range spanning between $~8$ and $2500 m$. We define the effective spatial range as the distance where the correlation equals 0.05 [18].

Three Markov Chain Monte Carlo (MCMC) chains were run for 20,000 iterations, each with the computationally demanding model requiring approximately 30 min to complete a single MCMC chain. We diagnosed convergence using the CODA library in the R Statistical and Computing environment by monitoring the mixing of chains and the Gelman–Rubin statistic [36]. Acceptable convergence was established within 10,000 iterations for all the models. The posterior inference was premised on a post-burn-in subsample of 15,000 iterations, that is, every third sample from the last 15,000 iterations of each chain. SVC and SVI models were fit using the spBayes R Statistical and Computing Library version 0.4.3. We, therefore, utilized the spBayes R statistical package for fitting all the predictive models for both Landsat-8 and Sentinel-2 SVC models.

2.6. Landsat OLI and Sentinel-2 MSI Imagery Derived Covariates

Landsat-8 has a revisit period of 16 days and offers nine spectral bands with a spatial resolution of 30 m for Bands 1 to 7 and 9 [6,37]. The panchromatic band, Band 8, has a spatial resolution of 15 m. On the other hand, Sentinel-2 has thirteen spectral bands where four bands are configured at 10 m spatial resolution, six bands at 20 m, and three bands at 60 m spatial resolution [37]. Common vegetation indices utilized in the estimation of biophysical variables of Absorbed photosynthetically active radiation (APAR), Leaf Area Index (LAI), and biomass are Normalised Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI) and Soil Adjusted Vegetation Index (SAVI) [5]. We utilized the same Vegetation Indices (VIs) in this study. Readily available and cost-effective new-generation sensors (Landsat-8 and Sentinel-2) with improved spectral and spatial resolution were therefore utilized in the modeling of C stock under the spatially varying coefficients assumption.

We obtained Landsat 8 imagery from the United States Geological Survey Earth Explorer (http://earthexplorer.usgs.gov) as data ready for analysis (ARDs) on the 20 of September 2020. All the datasets were riddled with cloud cover and shadow cloud cover limits set to smaller than 10%. We acquired Sentinel-2 cloud-free images on 20 September 2020 at the same time as the Landsat 8 OLI data collection covering the entire area coinciding with the dates Landsat-8 OLI was collected covering the domain of interest. The multispectral instrument is the main imaging instrument used for Sentinel-2, a push broom scanner that measures the terrestrial Top of the Atmosphere (TOA) reflectance in thirteen spectral bands, that is, 443 nm to 2190 nm. We derived Sentinel-2 spectral data as level-1C 12-bit automated TOA reflectance values. Orthorectification and pre-processing of the TOA-derived products were performed using the sen2r package of the R Statistical and Computing environment [38].

Soil Adjusted Vegetation Index ($SAVI$), Enhanced Vegetation Index ($EVI$), and Normalized Difference Vegetation Index ($NDVI$ were used as Landsat-8 and Sentinel-2 derived covariates for C stock prediction in the plantation forest of the eastern highlands of Zimbabwe. The literature is replete with studies that have utilized the aforementioned vegetation indices [3,39] as ABG predictors in biomass estimation. In this study, predictor variables supplied by both Landsat-8 and Sentinel-2 are used specifically for comparing SVC Bayesian hierarchical geostatistical models predicting C stock in a disturbed environment in Zimbabwe. Given information regarding C stock distribution at each location, we fit the model in Equation (1) $SVC$ with $p=2$. Hence, we have two processes, an intercept and one slope process relating to $NDVI$.

2.7. Model Fit and Prediction Accuracy Evaluation

We assessed the performance of the models using the commonly used Deviance Information Criterion (DIC) to categorize models in terms of how well they fit the data [40]. The sum of the Bayesian Deviance and the effective number of model parameters make up the DIC criterion. The Bayesian deviance measurement, which assesses model goodness of fit, and the effective number of model parameters which penalizes model complexity, are measured by D and pD, respectively. Attractive models have lower DIC values.

Predictive performance was evaluated through a $k-fold$ cross-validation technique. C stock was predicted from observations within each subset, given the estimated parameters from the remaining subsets. We employed the Root Mean Squared Prediction Error (RMSPE) as a metric from the R Statistical and Computing environment to calculate the sampled C stock data values and the accompanying median of the posterior predictive distribution (PPD). Models with lower RMSPE signify more accurate C stock predictions.

3. Results

3.1. Multispectral Remote Sensing C Stock Derived Predictors

Employed predictors in C stock prediction using a Bayesian hierarchical framework with spatially varying coefficients showed NDVI as a significant predictor of C stock as illustrated in

Table 1. Landsat-8 and Sentinel-2 derived predictors of C stock.

Parameter	Lansat-8 OLI C Stock Model				Sentinel-2 MSI C Stock Model
	Mean	s.d	2.5%	97.5%	Mean	s.d	2.5%	97.5%
$\mathit{I}\mathit{n}\mathit{t}\mathit{e}\mathit{r}\mathit{c}\mathit{e}\mathit{p}\mathit{t}$	−2.65	1.02	−4.76	−0.83	−2.40	0.31	−3.01	−1.80
$\mathrm{NDVI}$	2.47	0.98	0.78	4.62	4.90	0.23	4.52	5.38
$\mathrm{SAVI}$	−0.57	0.67	−2.04	0.63	−0.55	0.34	−1.25	0.17
$\mathrm{EVI}$	−0.62	0.55	−1.57	0.50	−0.002	0.096	−0.17	0.24
${\sigma }_w^2$	1.25	0.26	0.72	1.72	0.075	0.016	0.046	0.092
${\sigma }_{\epsilon }^2$	0.35	0.17	0.074	0.54	0.0043	0.0030	0.0005	0.011
$\varphi$	0.0016	0.0003	0.0012	0.0018	0.0015	0.0002	0.0014	0.0014

. 95% credible intervals of SAVI and EVI contained zero and hence, rendering them as insignificant predictors of C stock. These predictors were therefore excluded in the final prediction and mapping of C stock distribution.

3.2. Candidate Models and Parameter Estimates

Model posterior estimates of the regression coefficients for the Landasat-8 and Sentinel-2-based C stock models are illustrated in

Table 2. Landsat-8-based SVC model median parameter estimates alongside their 95% credible intervals.

		Non-Spatial	SVI	$SVC-NDVI$	SVC
Parameter C.I 50% (2.5%, 97.5%)	${\beta }_0$	−2.5 (−4.1, −0.8)	−2.7 (−4.6, −0.8)	−2.9 (−5.0, −0.9)	−3.5 (−5.3, −1.7)
	${\tilde{\beta }}_{NDVI}\left(s\right)$	5.4 (3.3, 7.5)	2.8 (0.8, 4.7)	3.3 (1.2, 5.1)	3.6 (1.4, 5.8)
	${\tau }^2$	1.5 (1.2, 2.0)	0.5 (0.3, 0.8)	0.12 (0.03, 0.40)	0.1 (0.02, 0.2)
	$\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{${\varphi }_0$}\right.$	-	0.0014 (0.0012, 0.0021)	0.003 (0.002, 0.003)	0.1 (0.07, 0.16)
	$\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{${\varphi }_{NDVI}$}\right.$	-	-	0.13 (0.06, 0.22)	0.03 (0.02, 0.05)
	${\sigma }_0^2$	-	1.0 (0.5, 1.5)	1.1 (0.8, 1.6)	0.24 (0.12, 0.48)
	${\sigma }_{NDVI}^2$	-	-	1.0 (0.2, 3.0)	0.45 (0.14, 1.43)
Fit Statistics	$D$	436	140	41.8	65.9
	$pD$	4.1	68.9	87.0	108.9
	$DIC$	207	105.6	−77.9	29.0
	$RMSPE$ (MgCha−1)	8	7.77	7.54	6.42

C.I means 95% Credible Interval.

and

Table 3. Sentinel-2-based SVC model median parameter estimates alongside 95% credible intervals.

		Non-Spatial	SVI	$SVC-NDVI$	$\mathbf{S}\mathbf{V}\mathbf{C}$
Parameter C.I 50% (2.5%, 97.5%)	${\beta }_0$	−2.5 (−4.1, −0.8)	−2.8 (−4.6, −0.9)	−2.9 (−4.9, −0.7)	−3.5 (−5.4, −1.6)
	${\tilde{\beta }}_{NDVI}\left(s\right)$	5.4 (3.3, 7.5)	3.0 (1.2, 4.9)	2.9 (0.7, 4.9)	3.7 (1.5, 5.7)
	${\tau }^2$	1.5 (1.2, 2.0)	0.5 (0.3, 0.9)	0.32 (0.15, 0.55)	0.2 (0.11, 0.42)
	$\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{${\varphi }_0$}\right.$	-	0.0015 (0.0014, 0.0016)	0.0016 (0.0015, 0.0018)	0.06 (0.018, 0.076)
	$\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{${\varphi }_{NDVI}$}\right.$	-	-	0.09 (0.036, 0.13)	0.0015 (0.0014, 0.0017)
	${\sigma }_0^2$	-	0.7 (0.4, 1.1)	1.1 (0.8, 1.9)	0.23 (0.16, 0.47)
	${\sigma }_{NDVI}^2$	-	-	0.7 (0.43, 1.1)	0.65 (0.41, 1.05)
Fit Statistics	$D$	436.6	159	93.3	17.9
	$pD$	4.1	58.0	89.8	106.1
	$DIC$	207.6	112.5	65.6	−177.3
	$RMSPE$ (MgCha−1)	7.85	7.69	7.46	6.23

C.I means 95% Credible Interval.

, respectively, for the non-spatial, SVI, and SVC models. Spatial autocorrelation is modeled in the residuals in the SVI, SVC, and all the $SVC-{\beta }_k$ variant models. This may entail differences in the posterior regression coefficient parameter estimates, depending on the C stock model parameter structure. Credible intervals (95% CI) for ${\beta }_0$ and ${\beta }_{NDVI}$ including zero for both non-spatial and SVI models would hint at a non-significant relationship between the C stock observations and predictor variables. However, we cannot apply the same reasoning and interpretation for the SVC models as the predictor-specific spatially varying coefficient maps of ${\beta }_{NDVI}+w_{NDVI}\left(s\right)$ should be considered and ascertain whether location-specific CIs include zero.

C.I means 95% Credible Interval.

C.I means 95% Credible Interval.

Posterior estimates for the Landsat-8 and Sentinel-2 based ${\beta }_{NDVI}+w_{NDVI}\left(s\right)$ coefficients are shown in Figure 3 and Figure 4, with a significant relationship between the outcome variable and $NDVI$ being evident over the entire study domain. In both cases, there is significant variability in the values of ${\beta }_{NDVI}+w_{NDVI}\left(s\right)$ in the studied region. This could imply more accessibility to plantation resources by communities settling in the plantation areas. The same trend is noticeable for the Sentinel-2-based C stock-based SVC model (Figure 5 and Figure 6) for both predictor coefficient values of ${\beta }_{NDVI}+w_{NDVI}\left(s\right)$ and for the same region as observed in the Landsat-8 based SVC model in Figure 3. Communities settled within certain areas of the plantation forest have more access to forest resources than in other areas, rendering low C stock density in these areas as demonstrated by the corresponding low $NDVI$ values in the same region for the Sentinel-2-based SVC model. Enhanced detail in the spatial resolution of Sentinel-2 based ${\beta }_{NDVI}+w_{NDVI}\left(s\right)$ vindicates variability in ${\beta }_{NDVI}+w_{NDVI}\left(s\right)$ coefficient values in the same region over those derived from generalized Landsat-8 multispectral data [41].

Estimates of ${\beta }_0+w_0\left(s\right)$ for the Landsat-8 and Sentinel-2-based SVI and SVC models are shown in Figure 3 and Figure 4, respectively. The ${\beta }_0+w_0\left(s\right)$ pattern for the SVI is the same for both Landsat-8 and Sentinel-2 spatially varying coefficient models whilst the ${\beta }_0+w_0\left(s\right)$ SVC pattern in Landsat-8 is the same for Landsat-8-based SVC (${\beta }_{NDVI}+w_{NDVI}\left(s\right)$) process. On the other hand, the ${\beta }_0+w_0\left(s\right)$ SVI pattern for Sentinel-2 based model is different for Sentinel-2 ${\beta }_{NDVI}+w_{NDVI}\left(s\right)$ SVC model. However, the same is not true for Sentinel-2-based SVI as the partitioning of $w_0$ into $w_0$ and $w_{NDVI}$ in Sentinel-2 SVC is detailed with enhanced spatial resolution compared to Landsat-8 based $w_{NDVI}$ SVC model.

Table 2 and Table 3 illustrate the posterior estimates of the uncorrelated residual variance, ${\tau }^2$. The uncorrelated residual variance is largest in the non-spatial model, whilst it is small for the SVC-variant models for both Landsat-8 and Sentinel-2-based SVC. The SVI and the SVC variant models incorporate a spatially varying correlated random effect, $w_0$ with variance ${\sigma }_0^2$. The SVC model variants further incorporate more spatially varying correlated random effects $w_{NDVI}$ and variance ${\sigma }_{NDVI}^2$. For both Landsat-8 and Sentinel-2 based models, the $w_0$ and $w_k$ explained much of the residual variability and hence, a reduced ${\tau }^2$. The implication is higher predictive accuracy for the SVC-variant models in both sensors, making the SVI and SVC more attractive over the error-independent models. This is supported by the goodness of fit diagnostics for the SVI and the SVC-variant models illustrated in Table 2 and Table 3 for the Landsat-8 and Sentinel-2 derived regression coefficients, respectively.

In comparison to the SVI model, the SVC model based on both remote sensing-derived covariates reduced the non-spatial residual spatial dependence by incorporating the space-varying impact of ${\beta }_{NDVI}$. Estimates of the spatial process parameters have a big difference. In particular, the spatial process parameters for the Landsat-8-based (Figure 3) SVC point estimates of ${\sigma }_0^2=0.24$ and the effective range of 30 m (i.e., $\approx -log(0.05/\varphi )=-log(0.05/0.1)$ denote a less variable and significantly shorter effective spatial range than the spatial process of the SVI model. The same pattern recurs in the Senitnel-2-based SVI and SVC models, where the effective spatial range reduces from $\approx 1800 m$ in the SVI to $\approx 100 m$ in the SVC model (Table 3). The non-negligible spatial process parameter estimates of ${\sigma }_{NDVI}^2$ and ${\varphi }_{NDVI}$ in the SVC model denote a potentially space-varying relationship between C stock and multispectral remote sensing derived covariates.

3.3. Landsat-8 and Sentinel-2 C Stock-Based Predictions

The entire SVC model based on Landsat-8 and Sentinel-2-based C stock models generated the lowest DIC, D, and RMSPE, as illustrated in Table 2 and Table 3, respectively. Landsat-8-based SVC 95% CIW is much shorter when compared to the non-spatial fitted model values. A 10% improvement in the RMSPE is seen in the Landsat-8-based model when moving from the error-independent to the SVC model. In the same vein, the Sentinel-2-based model had a 12% improvement from the SVI to the SVC model (Table 3). Predictions produced for both Landsat-8 and Sentinel-2-based C stock models mirrored the observed C stock data in the studied region. Observed C stock data values ranged from log (0.2–6) MgCha⁻¹ (Figure 1). Variability in C stock uncertainty is fairly constant across the studied domain for both Landsat-8 and Sentinel-2 SVC-based models.

The almost similar variability in the density of C stock predicted by both sensors could be attributed to the inadequacy of covariates in the modeling framework, which when the range of modeled covariates is broadened, could accurately depict the variability of C stock in these disturbed plantation forest ecosystems [5].

3.4. C Stock Model Prediction Assessment

Scatterplots observed against predicted C stock alongside the 95% intervals are illustrated in Figure 7. Slight improvement in the performance of Landsat-8 and Sentinel-2 C stock-based predictions can be validated through a visual inspection of the results. Estimated Root Mean Square Prediction Error (RMSPE) for Landsat-8 (6.42 Mg C ha⁻¹) and Sentinel-2 (6.23 Mg C ha⁻¹) based C stock prediction further reinforces the model prediction performance for the two sensors illustrated in Figure 7.

4. Discussion

Space varying coefficient models constructed from different but related new-generation multispectral remote sensing platforms were used to predict C stock in a managed plantation forest ecosystem in Zimbabwe. The RMSPE is marginally higher in Sentinel-2-based C stock SVC model than in the Landsat-8-based SVC counterpart. Our findings with regard to the performance of Landsat-8 and Sentinel-2 sensors are also congruent with the work of [12,37]. SVC models for both new-generation remote sensing-derived predictors showed preference over the error-independent models. However, estimates from the two data sources were marginally different from each other on the prediction of C stock, as illustrated in validation diagnostics [42]. The adaptable structure of the SVC permitted the residuals spatial variability to be apportioned between the random slope and the random intercept. This provided additional benefits not available in the SVI models from the two predictor sources. The SVC permitted reconnaissance of the C stock observations and the $NDVI$ predictor. Furthermore, the SVC models in both multispectral data sources had better representation of the processes thereby yielding C stock predictions with reduced variance.

Evaluation of the models utilizing predictors from both remote sensing data sources showed SVC to fit the data better than the SVI models. Kriged maps for C stock using Landsat-8 and Sentinel-2 data were not significantly different from each other, with the Landsat-8 SVC displaying a slightly wider 95% CI compared to the Sentinel-2-based SVC model. Again, this is partly because the study employed conventional bands (indices) that are calculated in a similar fashion in both Landsat-8 and Sentinel-2, and hence, the differences in prediction only emanated from the spatial resolutions. Explicit accommodation of residual spatial dependence through spatially correlated random effects gathered better predictions as the SVC models using the two data sources as regression coefficients borrowed additional information from neighboring C stock observations [43,44]. Precise estimates of covariance parameters are not easy to derive with small inventory sample sizes [45]. Consequently, we might anticipate some impact on the accuracy of the predictions when uncertainty in the covariance ensues to the posterior predictive distributions. Predictions had limited information to borrow from because of the sparsity of C stock inventory observations. This was further worsened by the reduction in the overall sample size through cross-validation.

Similar to previous studies predicting ABG and C stock, our findings establish Sentinel-2 as a better source of RS data for predicting C stock in disturbed environments [37] compared Sentinel-2 and Landsat-8 imagery for forest biomass prediction and showed Sentinel-2 outperforms Landsat-8 because of the enhanced spatial resolution in the former in comparison to Landsat-8. Most studies comparing Sentinel-2 and Landsat-8 for predicting AGB prefer Sentinel-2 over Landsat-8. This is further justified by the work of [13], who compared Worldview-3, Sentinel-2, and Landsat-8 for representing AGB in a forest environment in Thailand and demonstrated Worldview-3 and Sentinel-2 as better data sources and, therefore, predictors than Landsat-8 due to the red-edge and the improved spatial and spectral properties of Worldview-3 and Sentinel-2. Furthermore, [16] utilized LiDAR derived covariates from establishing the prediction performance of SVC models using forest inventory data in North America. The researchers established significant improvement in biomass prediction accuracy in the presence of residual spatial dependence deriving from the finer resolution LiDAR covariates. In most cases, the effectiveness of SVC models in these studies is usually strengthened by the solid non-stationary relationships between the response variable and the predictor variables influenced by unobserved ecological factors operating at broad geographical scales. Such ecological factors were also seen in the present study as NDVI was established to be a statistically significant predictors of C stock in the studied region. The biggest drivers to these factors are the presumed activities of forest disturbance due to human encroachment into plantation forests that subsequently impact the density and distribution of ABG biomass.

4.1. Limitations of the Study

Apart from the vegetation indices influencing the spatial distribution and density of AGB employed in this study, it is also known and acknowledged in the literature that climate and topographic variables play a part in the distribution of C stock. For example, [37,38] have shown elevation and aspect accounting for the bigger portion of the spatial variability of C stock in a mountainous region of Nepal. As such, the limitation of our study is the application of vegetation indices as predictors of C stock, and this may, therefore, not be representative of the general C stock dynamics in the studied region. We, therefore, recommend the integration of topo-climatic factors with new generation remote sensing-derived vegetation indices for future research in order to obtain a more accurate global overview of the C stock density and distribution in the studied region.

4.2. Conclusions

The study presented a hierarchical Bayesian geostatistical spatially varying coefficient model for determining the relationship between sampled C stock data and multispectral remote sensing derived predictors. There was a marginal improvement in model fit, and prediction accuracy in both Landsat-8 and Sentinel-2-based SVC models in comparison to the error-independent models. The SVC model permitted exploration of the observed C stock locations where the models performed well or poorly, which was missing in the SVI models. This provided an understanding of the performance of the multispectral remote sensing derived predictors for modeling C stock and hence, sets the foundation for the updating of the carbon forest plantation database for forest practitioners in the country and utilized as a monitoring tool in the long term. The Sentinel-2-based SVC model was preferred for prediction in the plantation forest ecosystems as its model provided tighter credible intervals compared to the Landsat-8-based C stock SVC model.

The small sample size of the data utilized in the present research enabled the modeling approach to be computationally feasible. When inventory plots are comprised of bigger sample sizes, the matrix operations of immense dimensions are needed for computing model parameter estimates of higher magnitude and may not be achievable through ordinary PCs. Our future work is therefore aimed at exploring algorithms for resolving the dimensionality curse when fitting spatially varying coefficient models. The problem of dimensionality can also get complicated if many predictors are involved in the SVC modeling framework. Resolving dimensionality issues is needed as forest C stock is typically modeled with predictors from many data sources, chief amongst them being topographic, bioclimatic, and anthropogenic variables.