Spatial variability of soil carbon across a hillslope restoration planting in New Zealand

Study area

The study was carried out near Glorit (174.4530E, 36.4921S; 174.4587E, 36.4950S), Auckland, New Zealand (Fig. 1). The 24-ha study site is classified as undulated medium hill area. The elevation ranges from 20 to 136 m above sea level with slopes reaching 20 degrees. The region where the site is located features a warm subtropical climate, with annual rainfall of 1457.5 mm and mean annual temperature of 14–16°C (Chappell 2014). The geology of the study site is late quaternary alluvium and colluvium, which consists of unconsolidated to poorly consolidated sand, gravel, mud, and peat. This overlies the early Miocene Waitemata group, which consists of interbedded sandstone and mudstone (Ballance and McCarthy 1975). The soil types of this area are classified as mottled and typic yellow ultic soils (Hewitt 2010).

Fig. 1.

Map showing the location of the study area (outlined in red), Kaipara Harbour, Rodney District, North Island, New Zealand.

During the late 19th and early 20th century, most of the dense forest and bush was cleared and converted into pasture (Beever 1981). Some areas of remnant bush remained and there are patches of naturally regenerating bush across the hillslope (Noble and Harmsworth 1990). In 2008 restoration work began and native species were planted annually until 2016.

Sampling design

The 121-sample points were chosen to adhere to budget and time constraints. A combination of systematic and judgemental sampling was used to generate sample points across the study area (Fig. 1). Systematic sampling based on a uniform grid was initially used, as it is preferable for spatial analysis. However, judgemental sampling, refined from a stratified random method was incorporated. This adjustment was necessary because certain vegetation groups were not adequately covered with sample points using systematic sampling alone (Lark and Marchant 2018), or roads hindered systematic coverage of the study site. The study area was divided into seven different vegetation groups (Fig. 1). The vegetation groups were pasture (sample points = 3); plantings between 2008 and 2010 (sample points = 18), 2011 (sample points = 13), 2012–2013 (sample points = 17), 2014–2015 (sample points = 35), 2016 (sample points = 6) and a reference group of bush (sample points = 29) that included the remnant bush and naturally regenerating bush present in the study area.

Vegetation survey and soil sampling

A vegetation survey was conducted at each sample point using the Point-Centred Quarter (PCQ) method (Mitchell 2010). Organic layer (forest floor) samples were taken using a PVC ring (diameter, 10 cm; length, 8 cm) at 121 locations. Mineral soil samples (one sample for bulk density and one sample for soil physical and chemical analysis) were collected using a soil auger (diameter, 3.5 cm). At 0–10 cm depth, samples were taken at all 121 locations, at 10–20 cm, samples were taken at 50 locations, and at 20–30 cm, samples were taken at 25 locations. Field measurements and sampling were conducted between October and December 2019.

Laboratory analysis

Organic layer samples were oven dried at 60°C for 1 week and then weighed to determine the organic layer biomass. Mineral soil samples for bulk density analysis were dried at 105°C for 1 week. The soil weight (oven dry) was divided by the moist volume to obtain bulk density (Blakemore et al. 1987). Mineral soil samples for physical and chemical analyses were dried at 40°C for 1 week (McKenzie 2000). Roots and organic matter were removed, and the soil was sieved through a 2-mm sieve. For particle size analysis 20 mL of 0.55% Calgon deflocculant solution was added to 10 g of soil. The soil–Calgon mix was analysed using a laser particle size analyser (Mastersizer 3000, Malvern Panalytical, Worcestershire, UK). For pH, 25 mL of deionised water was added to 10 g mineral soil and measured using a pH electrode (sensION + PH3, HACH, Loveland, CO, USA).

Accurate and precise determination of soil carbon concentration can be obtained by automated, high-temperature dry combustion methods (Chatterjee et al. 2009). Due to financial constraints and instrument availability, we were only able to analyse a subset of mineral soil samples (n = 66) by high-temperature dry combustion. The oven-dried (<2 mm) samples were ground into powder using a centrifugal mill (ZM 300, Retsch GmbH, Haan, Germany). Total carbon concentration was determined by weighing 10 mg of ground material into tin cups and combusting the sample at 950°C using an elemental analyser (Vario El Cube Elementar, Langenselbold, Germany). Of these, 10% of samples were measured twice. The variability was <1%.

Loss on ignition (LOI) is a cheaper method to estimate organic matter content from which the total carbon concentration can be calculated using a conversion factor. When LOI and total carbon concentrations are both measured, regression models are used to convert LOI to soil carbon concentration (Jensen et al. 2018). We determined organic matter in all mineral soil samples by LOI following the method described by Blakemore et al. (1987). Five grams of oven-dry (40°C) soil was placed into previously ignited and weighed porcelain crucibles, dried at 105°C for 24 h, cooled in a desiccator and weighed again. Finally, the crucibles were ignited at 500°C for 4 h in a muffle furnace (Elecfurn, The Electric Furnace Co Ltd, Auckland, New Zealand). After ignition, the crucibles were cooled in a desiccator and weighed. The organic matter content was calculated as the difference between the oven-dry (105°C) weight before and after ignition.

A linear regression function was established between organic matter content and total carbon concentration (r² = 0.88; P < 0.001). This function was used to estimate total carbon concentration for the remaining mineral soil samples. The pH values (median = 5.3) suggest an insignificant amount of inorganic (carbonate) carbon. Thus, total carbon concentration is a good reflection of soil organic carbon in this study area. Soil carbon stocks (kg m⁻²) were calculated using the estimated total carbon concentration (g kg⁻¹ soil), bulk density (kg m⁻³), and sampling depth (m).

Landscape variables

To calculate elevation, slope, and aspect, we used an 8-m Digital Elevation Model (LINZ 2023a). The road and river centrelines, essential for distance calculations concerning paths and water bodies, were sourced from the NZ Road/Rivers Centrelines dataset (Topo, 1:250 k) (LINZ 2023b). To ensure precision in distance measurements, we generated Euclidean distance surfaces. The Euclidean distance surface illustrates the distance from each raster cell to a source or a set of sources (rivers, roads) based on straight-line distance.

Statistical analysis

Shapiro-Wilk’s test was used to test if the data was normally distributed (Shapiro and Wilk 1965). As some soil characteristics were non-normally distributed, a Kruskal–Wallis test was used to identify if there were significant differences in soil characteristics between vegetation groups at P < 0.05 (Kruskal and Wallis 1952). A post hoc Wilcoxon test was then run to identify which vegetation groups were significantly different. Statistical tests were carried out using R (R Core Team 2021).

Spatial analysis

Spatial interpolation and MGWR methods were used to predict the spatial distribution of soil carbon stocks (0–10 cm, 10–20 cm, and 20–30 cm) and to investigate the relationship between soil and vegetation characteristics and topographic features for the top layer of soil (0–10 cm). The outline of the analysis is in Fig. 2.

Fig. 2.

Flow diagram outlining the spatial analysis.

To understand the spatial distribution of soil carbon stocks, ordinary kriging was selected over other univariate geostatistical methods as this method accounts for local variation in a studied phenomenon and so predicts spatial trends in values of an attribute (Li and Heap 2011, 2014). The accuracy of this method is limited, because it uses a stationary mean within the local search window. To execute ordinary kriging, the spatial dependence of observations was specified through a variogram. An experimental variogram is determined directly from measurements of the characteristic, and the prediction of unsampled locations is then generated by fitting a theoretical variogram to the experimental variogram (Volfová and Šmejkal 2012). Theoretical variograms provide continuity so variance can be obtained at any distance, and follow various model forms (i.e. spherical, exponential, steady/stable). A variogram contains a sill, nugget, and range. The sill is the upper value of the variogram, the nugget is the value of the variogram where the distance is zero, and the range is the maximum distance at which two observations are dependent (Volfová and Šmejkal 2012). To assess how accurate the predicted values at unknown locations were, leave-one-out cross-validation was performed where one data point was held back from analysis and predicted using the data from all other locations. The predicted and the original values are then compared. This is then repeated with all data points to inform the quality of the model. One operation of cross-validation was calculating the root mean square error (RMSE). Spatial interpolation was carried out using ArcMap ver. 10.8 and variograms were generated using ‘gstat’ and ‘automap’ packages in R (Environmental Systems Research Institute 2020; R Core Team 2021; Himestra and Skoin 2023; Pebesma and Graeler 2023).

To explore the spatial variability in the relationships between soil carbon stocks and potential influencing factors, we employed a multiscale geographically weighted regression (MGWR) analysis.

A standard geographically weighted regression (GWR) represents a spatial statistical approach that surpasses conventional global regression models, which assume constant relationships across space. It achieves this by acknowledging the spatial heterogeneity in spatial processes, which delineates how these processes vary according to spatial context (Fotheringham et al. 1998; Oshan et al. 2019). In GWR, local linear models are established at multiple locations using data from the neighbourhoods surrounding these points within a specified spatial bandwidth. Multiscale geographically weighted regression considers the spatial coordinates of the observed phenomena and predicts parameter estimates based on weighted distances between these locations and a given regression point (Fotheringham et al. 2017). Moreover, MGWR extends beyond GWR by accounting for the spatial influence on the relationship between independent and dependent variables. It allows for the modelling of both local and global relationships for all variables within the model. In practical terms, MGWR computes separate regression models for every location within the dataset, and it calculates a corresponding bandwidth for each variable. This model can be formally expressed as:

where y_i is a value of the dependent variable at i location, β₀ is an intercept, β_bwj is the bandwidth to calibrate the ith relationship, (u_i, v_i) is the spatial location of the ith observation, and ε is the random error term.

To enhance our understanding of the relationships between soil carbon stocks and various characteristics encompassing vegetation, soil, and site attributes, we employed a MGWR model. We systematically selected 13 exploratory variables, guided by the existing literature and data availability. Subsequently, we conducted a model selection procedure based on the Akaike Information Criterion (AIC) (Fotheringham et al. 2002). This rigorous process uses AIC values to identify suitable variables based on stepwise forward selection process. Observations revealed that after incorporating a total of 11 variables, the AIC values ceased decreasing with the addition of further variables. Consequently, the last two aspect-related variables were excluded from the final model. As the dependent variable we used soil carbon stocks from 121 sampling locations.

The MGWR was performed using a Python library developed by Oshan et al. (2019) and the results were visualised in ArcMap ver. 10.8.

Landscape, vegetation, and soil characteristics

Stem density (0.1–2.9 trees m⁻²), stand basal area (0.2–645 m² ha⁻¹) and organic layer biomass (1–104 g m⁻²) varied considerably across the study area (Table 1).

The soil was sand and silt dominated and classified as loamy sand, sandy loam, and loamy silt. Median bulk density ranged between 1 and 1.3 g cm⁻³ and increased slightly with soil depth. Median pH at 0–10 cm and 10–20 cm was 5.3 and slightly lower at 20–30 cm (Table 1).

Soil carbon stocks at 0–10 cm varied between 1.9 and 7.1 kg m⁻² (Fig. 3). Soil carbon stocks in the pasture was significantly higher than in restoration plantings and bush. Soil carbon stocks in the 2014–2015 plantings were higher than 2011. No difference in soil carbon stocks was found between the other vegetation groups. Carbon stocks at 10–20 cm depth ranged from 1.8 to 5.4 kg m⁻², and in 20–30 cm depth from 2.3 to 4.4 kg m⁻². There were no significant differences in soil carbon stocks between vegetation groups at 10–20 cm and 20–30 cm.

Fig. 3.

Soil carbon stocks (0–10 cm) across the vegetation groups. The horizontal line in the box indicates the median; the boundaries of the box represent the 25th and 75th percentile, and the whiskers show the highest and lowest values.

Soil carbon stocks found in the bush and pasture at Glorit were within the range found in New Zealand indigenous shrublands (3.2–4.6 kg m⁻²; Hedley et al. 2013) and pastures (3.7–6.7 kg m⁻²; Scott et al. 1999; Sparling and Schipper 2004). Soil carbon stocks varied considerably within a given layer but tended to decrease with depth (Fig. 4). A decrease in soil carbon stocks with depth has been reported in other studies and is partly explained by the decline in litter and root derived organic matter input with depth (Lawrence et al. 2015; Henneron et al. 2022).

Fig. 4.

Ordinary kriging of soil carbon stocks (kg m⁻²) at different depths across the study area. (a) 0–10 cm, (b) 10–20 cm, and (c) 20–30 cm. Note that we did not map carbon stocks from pasture-related samples as we only had three sampling locations.

Soil carbon stocks (0–10 cm) in the restoration plantings were significantly lower compared to the pasture (Fig. 3). This is in line with previous findings (Guo and Gifford 2002; Guo et al. 2007) reporting a decrease in soil organic carbon following afforestation. The young age of the tree plantings may partly explain the lower carbon stocks in the restoration plantings (Guo and Gifford 2002). Another factor may be the initial carbon stock. An increase in soil carbon following afforestation often occurs in soils low in carbon, while a decrease has been found in carbon rich soils (Hong et al. 2020).

Ordinary kriging

We computed and fitted a variogram to assess the distribution of soil carbon stocks across various vegetation groups. Pasture samples were excluded from the analysis due to their limited number (three in total) and their spatial separation from the primary study site. The experimental variograms of soil carbon stocks (0–10 cm) had relatively weak trends, which can happen in complex terrain with numerous variables at play. The theoretical models fit the experimental variograms with the following parameters (model choice, steady; sill, 0.45; Nugget, 0.43; Range, 66.11).

The predicted surface of soil carbon stocks at 0–10 cm showed a small area of low values (2.5–2.8 kg m⁻²) at the centre of the study area (Fig. 4a). Higher values (4.4–6.2 kg m⁻²) were predicted at the south-western edge and the southern tip at the eastern end of the study area. The predicted surface of soil carbon stocks at 10–20 cm showed lower values (2.9–3.3 kg m⁻²) in the centre, and higher values (3.4–4.4 kg m⁻²) in the north-west of the study area (Fig. 4b). The prediction surface for this layer was smoother than that at 0–10 cm depth. At 20–30 cm low values (2.3–2.8 kg m⁻²) were again found in the centre of the study area and high values (4.0–4.4 kg m⁻²) were predicted in the north-western section of the study area (Fig. 4c).

Multiscale geographically weighted regression

Results derived from the calibrated MGWR model are in Table 2 and Fig. 5. Table 2 shows the final list of variables with ranges of parameter estimates as well as model’s diagnostics.

Fig. 5.

Multiscale geographically weighted regression (MGWR) local parameter estimates for (a) elevation, (b) slope, (c) stem density, (d) bulk density, and (e) pH.

Fig. 5 shows five maps of local parameter estimates. Only local parameters that were significant were mapped. The global estimate demonstrated a statistically significant negative relationship (−0.426), signifying an inverse relationship between elevation and soil carbon stocks (Fig. 5a). This suggests that carbon stocks decreased with elevation, which can be attributed to slightly decreasing temperatures, resulting in reduced rates of organic matter decomposition (Saeed et al. 2019). This trend was particularly strong in the western site of the study area, characterised by lower elevation and situated at the base of the hill. Lower elevation in this study area coincided with comparatively low bulk density values and low basal area (see Supplementary Figs S1–S4), which may partly explain lower soil organic carbon stocks. The inverse distance interpolation (IDW) was chosen to map chemical, physical, and environmental parameters shown in supplementary materials. This method calculates value in a grid point while considering neighbouring locations within a user-defined search radius. IDW is commonly employed in soil investigations due to its ease of implementation (Li and Heap 2011, 2014; Zhang et al. 2020).

Fig. 5b illustrates the local parameter estimates concerning the association between soil carbon stocks and slope, accounting for variations in other variables within the model. The global parameter estimate for this relationship was −0.258, suggesting an overall negative influence of slope on soil carbon stocks. However, the results from the MGWR indicate spatial variability in the strength of these effects. Across most of the western and central portions of the study area, the parameters were statistically insignificant. In contrast, the eastern part showed a significant negative relationship, indicating that increasing slopes in this region corresponded to lower soil carbon stocks. A negative relationship between slope and soil carbon has been reported in previous studies and has been explained by soil erosion and water runoff washing organic material off steep slopes to lower elevations (Gregorich et al. 1998; Zhong et al. 2017).

Fig. 5c depicts the local parameter estimates relating to a vegetation characteristic, stem density. The global relationship was positive (0.247), suggesting that higher stem density has a positive influence on soil carbon stocks. These parameter estimates indicate spatial variation, with significance observed exclusively in the western area of the study site with older (2008–2010) and more recent plantings (2024–2025). Some of the restoration plantings contained patches of higher stem density (Figs S1–S4). This is similar to the results from Na et al. (2021) where they found that higher tree density resulted in higher soil carbon stocks.

Fig. 5d indicates the study area where bulk density (g cm⁻³) had a significant positive effect on soil carbon stocks. The global estimate was significantly positive (0.709). The significantly positive relationship was found across the whole study area, with a stronger effect in the western area of the study site where bulk density tended to be lower compared to the remaining area (Figs S1–S4). Given that soil carbon stocks are estimated using soil carbon concentration and bulk density, bulk density was the dominant driver in areas where soil carbon concentration was fairly homogenous as in the western part of the study area (Figs S1–S4).

Fig. 5e delineates the local relationships between pH and soil carbon stocks. Globally, this relationship was not significant (0.142). However, the results from MGWR indicate local variations, with a significant positive relationship between pH and soil carbon stocks observed in the central part of the study site. Some studies have reported a positive relationship between pH and soil carbon content (Kuśmierz et al. 2023). Higher pH have been shown to contribute to enhanced plant growth resulting in higher biomass and litter input (Malik et al. 2018).

Further, we used MGWR to predict levels of soil carbon stocks in the study area (Fig. 6). We used IDW to map the results as a continuous surface.

Fig. 6.

Multiscale geographically weighted regression (MGWR)-derived soil carbon stocks (0–10 cm). These predicted soil carbon stocks values were interpolated using IDW interpolation to produce this continuous surface.

The MGWR-derived carbon stocks (Fig. 6) exhibit similar patterns to those derived from ordinary kriging (Fig. 4a). Predictions from MGWR reveal more detailed patterns with sharper boundaries for minimums and maximums. The fact that the predictions extend beyond the original values suggests that landscape and vegetation characteristics influence soil carbon stock levels. Typically, in flat areas, ordinary kriging may predict better due to accounting for spatial auto-correlation. However, in our case study, with a significant effect of elevation and slope, it is reasonable to assume that MGWR captures some variations more effectively (Wang et al. 2012).

In order to understand potential disparities in soil carbon stock values generated by the two selected methods, we have created a map illustrating the variations in soil carbon stock values stemming from distinct MGWR-derived predictions and kriging-based estimates (Fig. 7). In regions with fewer data points, there is a notable trend of larger negative differences, indicating that kriging tended to overestimate values while GWR underestimated them. Conversely, areas with higher bulk density exhibit a tendency toward larger positive differences, with GWR surpassing kriging estimates.

Fig. 7.

Map illustrating the variations in soil carbon stock values stemming from distinct MGWR-derived predictions and kriging-based estimates. Negative areas correspond to places where MGWR model predicted lower values of carbon stocks than kriging.