Assessing and Modelling Soil Organic Carbon Dynamics in a Mining Spill Restoration Area from 1998

1. Introduction

Soil organic carbon (SOC) is a fundamental component of terrestrial ecosystems, playing a crucial role in soil health, agriculture, the global carbon cycle, and climate change. It has been estimated that SOC accounts for around 25% of the potential of natural climate solutions, although approximately 60% of this potential corresponds to the restoration of depleted stocks [1]. Multiple benefits arise from soil carbon sequestration, including (1) the formation, protection, and remediation of soils and sediments, (2) regulation of risks and extreme events, (3) climate regulation, and (4) the creation and maintenance of habitats; all of which are aligned with the United Nations Sustainable Development Goals (SDGs) [2].

Despite this potential, soils cannot accumulate carbon indefinitely. In soils functioning as carbon sinks, SOC increases only until a new equilibrium is reached [3]. For instance, grassland ecosystems cannot serve indefinitely as carbon sinks, since their accumulation tends to stabilise over time under constant management and environmental conditions [4]. Acknowledging this limitation, the European Union Soil Strategy for 2030 explicitly emphasises the enhancement of SOC to promote soil health, biodiversity, and climate neutrality by 2050 [5]. Initiatives such as “4 per 1000” reinforce the message that even modest annual increases in global SOC can substantially offset emissions. However, achieving these political targets requires robust monitoring, reporting, and verification (MRV) systems for SOC changes, which cannot be achieved solely through conventional field sampling.

In this regard, remote sensing has emerged as a powerful tool for indirect SOC monitoring. Traditional methods are labour-intensive and spatially limited, whereas satellite data provide frequent and consistent coverage over large areas. Remote sensing enables the estimation of variables such as land cover, vegetation productivity, and even soil properties like clay content, all of which influence SOC dynamics. Recent studies, for example, have employed satellite-derived proxies of vegetation productivity and soil texture as inputs in SOC models, enabling broader spatial assessments [2,6]

Complementing observational tools, process-based models such as RothC are critical for simulating SOC dynamics over the long term. RothC (the Rothamsted Carbon Model) is designed to estimate SOC turnover in non-waterlogged soils, considering climate, soil texture, and organic inputs on a monthly basis. The model divides SOC into conceptual pools with different decomposition rates, allowing predictions of how carbon stocks will evolve under contrasting scenarios [7]. Recent applications have further extended RothC towards spatially explicit and hybrid modelling frameworks, combining process-based SOC simulations with satellite-derived indicators, gridded climate data and digital soil information to support regional carbon accounting and carbon farming assessments [8,9,10]. This type of input-substitution approach is particularly relevant because recent RothC applications have shown that model outputs can be strongly affected by the estimation of carbon inputs, evapotranspiration, initial SOC stocks and soil-property layers, especially when spatial datasets are used instead of locally measured information [8,10,11]

This study investigates SOC dynamics in the Guadiamar Green Corridor in southwest Spain, an area restored after a major mining disaster in 1998. Following the toxic sludge spill, large-scale remediation and land-use conversion efforts transformed the contaminated lands into a protected ecological corridor with native vegetation [12]. This unique landscape provides an ideal case to examine SOC evolution during ecological restoration processes, where huge clayley and organic amendments were used in the most polluted plots.

We analyse SOC using field campaigns conducted in 2000, 2014 and 2024 to assess temporal changes over a 25-year period. Simulations with the RothC model were run under different boundary conditions, comparing inputs derived from in situ measurements with those obtained through remote sensing. By assessing the consistency of model outputs with observed SOC values, we seek to determine whether satellite-derived parameters can reliably substitute field measurements in SOC modelling.

Beyond the plot-scale simulations, this study also explores SOC dynamics at landscape scale through a spatially explicit implementation of RothC driven by gridded environmental datasets. In this context, global baseline soil products such as SoilGrids were used to initialise SOC and clay conditions across the entire Guadiamar basin, allowing the evaluation of regional SOC dynamics under homogeneous spatial forcing. This landscape-scale approach provides an opportunity to assess how differences between locally measured baseline conditions and globally modelled soil inventories influence long-term SOC simulations, uncertainty propagation, and carbon balance estimates. Such comparisons are particularly relevant in restored or heavily disturbed ecosystems, where historical management interventions and anthropogenic soil reconstruction may not be adequately captured by global digital soil mapping products.

Within the broader framework of carbon accounting and climate change mitigation, monitoring, reporting and verification (MRV) systems are essential to ensure the credibility and robustness of national and international initiatives aimed at enhancing SOC [2]. These systems increasingly combine field observations, remote sensing, and process-based models to provide scalable assessments across multiple spatial and temporal scales. In this regard, the systematic replacement of field-derived variables by remote sensing proxies or open-access datasets offers a promising strategy to evaluate uncertainty and develop cost-effective Tier 3 MRV methodologies in data-scarce regions [5,13].

The main objectives of this study were to model the carbon dynamics in the RothC model to estimate changes of SOC despite of lack of in a Mediterranean climate, more specifically in a restored area as the Guadiamar Green Corridor. By using RothC, we aimed at including a strategy based in boundary conditions to evaluate the feasibility of replacing in situ measurements with estimates derived from remote sensing in SOC simulations using RothC mode. We aim not only to validate its accuracy and reliability but also to explore its potential to expand SOC monitoring in plot and landscape scales where field data are scarce, absent, or problematic [2,14,15].

2. Materials and Methods

2.1. Experimental Plots and Data Availability

The Guadiamar Green Corridor was established following the accidental mining tailings spill in Aznalcóllar on 25 April 1998, as an ecological recovery measure. Acid waters and sludge containing heavy metals from the tailings pond were discharged into the Agrio River and subsequently into the Guadiamar River. This event also caused both rivers to overflow, leading to widespread contamination of their floodplains [16]. The affected area covered 4,286 ha, of which 2,557 ha were agricultural lands [16].

The deposition of tailings on the soil not only introduced severe contamination from mining residues but also removed parts of the soil profile, reduced fertility, and altered its physical structure [17]. To remediate the affected soils, amendments were applied as an assisted natural remediation technique aimed at accelerating ecosystem recovery [18]. Inorganic amendments were used for arsenic stabilisation, while organic amendments were applied to restore soil fertility. Complementary techniques of phytoremediation and phytoextraction were also implemented [17].

The area is characterised by a Mediterranean climate, with mild, rainy winters (mean annual precipitation approximately 500 mm) and hot, dry summers. The mean annual daily temperature is around 17 °C, with a maximum of 33.5 °C in July and a minimum of 5.2 °C in January [19]. In the Guadiamar floodplain, soils generally present a neutral to slightly alkaline pH. However, in some terraces on the right bank, where gravel deposits predominate, lower pH values have been recorded. Soil texture ranges from sandy loam to clay loam [20].

Six long-term monitoring plots were considered in this study (Figure 1). These plots were established between 1998 and 2000 as part of the post-spill monitoring network and were considered landscape units representative of the same land-use or vegetation-cover type within the Guadiamar Green Corridor. The selected plots represented the main vegetation-cover types currently present in the study area: riparian forest (GUAD1, GUAD3 and GUAD6), shrubland (GUAD4 and GUAD5), and grassland (GUAD2). In this manuscript, the term “plot” refers to each of these six fixed landscape-monitoring units, while “soil replicates” refers to the individual soil samples collected within each plot during each sampling campaign.

In the first sampling campaign (2000), five soil replicates were collected from each monitoring plot at two depths (0–15 cm and 15–30 cm), resulting in a total of 10 samples per plot. Replicates were spaced approximately 15–30 m apart, with the aim of capturing the internal spatial variability of each plot while maintaining consistency with the historical monitoring design. In the 2014 campaign, three replicates per plot were again collected at the same two depths. The 2024 sampling campaign was designed to reproduce, as closely as possible, the previous sampling locations using the georeferenced points from the 2000 and 2014 campaigns. In contrast to the previous campaigns, the 2024 sampling focused exclusively on the upper 0–20 cm soil layer, since the main objective was to assess recent SOC dynamics and budget limitations prevented the collection of a larger number of samples and deeper soil profiles.The monitoring design was therefore based on long-term georeferenced landscape units rather than on a fully balanced experimental design. Consequently, the grassland class, represented by a single plot, should be interpreted as a site-specific long-term monitoring case within the restored corridor, rather than as a statistically replicated representation of all grassland areas.

SOC concentrations were determined using the Walkley–Black wet oxidation method [21]. Due to the application of organic–clay amendments and the dynamic nature of the soils, substantial differences were found in SOC contents, including abnormally high values compared with other Mediterranean Fluvisols [22]. Clay content was determined by particle-size analysis using the Boyoucos hydrometer method on field samples. For each study area, mean clay values were considered, since the application of amendments (such as red clays) altered soil granulometry in some zones of the study area, producing marked variability over short distances, like that observed in SOC concentrations.

2.2. RothC Configuration

In this study, we employed the RothC model (Rothamsted Carbon Model, version 26.3), developed by Rothamsted Research and described by Coleman and Jenkinson. The choice of RothC was based on several reasons: it requires a relatively small set of input data, it has been validated under a wide range of climatic and edaphic conditions, and it is widely used internationally for estimating SOC stocks, both in its original configuration and in adapted versions. In this case, the model was implemented in RStudio using the SoilR package [9,23,24].

The model structures SOC into five compartments: four active pools: decomposable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO), and humified organic matter (HUM), and one inert pool (IOM) [7,25]. Organic matter inputs to the soil are partitioned between DPM and RPM according to residue quality and land use. Carbon in these pools flows into BIO and HUM, with part being released as CO₂ through microbial decomposition processes. BIO and HUM also exchange carbon and release CO₂. Each active pool is characterised by a specific decomposition rate: DPM (0.165 years), RPM (2.31 years), BIO (1.69 years), and HUM (49.5 years). The total CO₂ flux resulting from these processes can be compared with the monthly heterotrophic respiration of soils, providing an estimate of soil carbon turnover under varying climatic and management conditions [7,26].

RothC requires a set of essential parameters to initialise soil conditions and define the factors regulating carbon dynamics. Among these, the initial SOC content serves as a reference for projecting its temporal evolution, while monthly carbon inputs—derived from plant residues (both above- and below-ground) or organic amendments such as manure or compost—represent the main source of decomposable organic matter [27]. Monthly climate variables, including precipitation, mean temperature, and potential evapotranspiration, further regulate microbial activity and thus the rate of carbon mineralisation [7,26]. Finally, clay content, expressed as weight percentage, plays a decisive role in SOC stabilisation, since higher fractions enhance its physical and chemical protection, thereby prolonging its persistence in the soil [29].

In the context of this study, several of these variables such as SOC, clay content, and climatic factors were obtained both from direct field measurements and from remote sensing products and statistical modelling. This combination enabled an evaluation, under the concept of boundary conditions, of the extent to which parameters estimated from remote sensing can reliably substitute for in situ measurements, providing a framework applicable to contexts where field data availability is limited. For the execution of the model, eight different scenarios were defined, referred to as boundary conditions (BC), and labelled from BC0 to BC7. These scenarios were linked to the five key parameters considered in RothC: soil organic carbon, clay, precipitation, temperature, and evapotranspiration.

2.3. Input Data

To run RothC we must have a set of certain variables: clay (%), soil organic carbon (t C ha-1), initial C stocks of the different pools (t C ha-1): DPM, RPM, BIO, HUM, IOM, average monthly precipitation (mm), average monthly mean air temperature (ºC), monthly open pan evaporation (mm)/evapotranspiration (mm) , yearly carbon inputs from plant residue (aboveground + roots + rhizodeposition),(t C ha-1) and other management variables showed in the Table 1.

The initial content of soil organic carbon (SOC) and clay was determined from sampling campaigns conducted in the study area between 1998 and 2000. To establish the baseline scenario, data from the year 2000 were used (Table 2) , collected only a few years after the restoration of the Guadiamar Green Corridor. The 2024 estimates are derived from new field campaigns applying consistent criteria regarding sampling depth, sample treatment, and analytical methodology to ensure temporal comparability. In this case, the SOC value employed corresponded to the mean of three replicates from the same plot. This approach was adopted due to the high variability in carbon content observed in this area because of organic amendment applications. Thus, for each study zone a single value was used, derived as the average of all collected samples (Table 2).

Annual carbon inputs represent the continuous incorporation of organic matter into the soil, either in the form of plant residues (both aboveground biomass and roots) or organic amendments such as manure or compost. As direct records were not available for the entire study period, these inputs were estimated by combining land cover data, residue production coefficients reported in the literature for each vegetation type, and assumptions about the proportion of biomass returned to the soil after the production cycle. To obtain realistic estimates of the annual carbon inputs required by the RothC model, an empirical approach was adopted that integrates vegetation dynamics observed through remote sensing with ecophysiological parameters drawn from the literature.

For each land use type (grassland, shrubland, and riparian forest), mean annual NDVI values were calculated from Landsat imagery (2000–2024), as shown in Figure 2. These values were then converted into the fraction of absorbed photosynthetically active radiation (fPAR) using the widely applied linear relationship fPAR = 1.24·NDVI – 0.168 [30]. Absorbed photosynthetically active radiation (APAR) was derived by multiplying fPAR by local estimates of incident annual photosynthetically active radiation (≈2000 MJ m⁻² yr⁻¹ in southern Spain, Running et al., 2000). Net primary productivity (NPP) was then estimated as the product of APAR and light use efficiency (ε), using ecosystem-specific values reported in the literature: 0.85 g C MJ⁻¹ for grasslands [32], 0.65 g C MJ⁻¹ for Mediterranean shrublands [31,33], and 0.95 g C MJ⁻¹ for riparian forests [34].

To approximate the fraction of NPP returning to the soil via roots, litterfall, and belowground detritus, specific allocation coefficients (fsoil) were applied: 0.40 for grasslands, reflecting their high root turnover and herbaceous inputs; 0.25–0.30 for shrublands, consistent with their more woody and recalcitrant biomass; and 0.15–0.20 for riparian forests, where a large share of NPP accumulates in persistent tissues [35]. These coefficients were further adjusted to reflect Mediterranean eco-physiological constraints such as summer drought and seasonal growth pulses. Finally, the fractions of rapidly decomposable and resistant materials were parameterised in RothC using the recommended DPM/RPM ratios for each land use (DR = 1.44 for improved grasslands, 0.67 for shrublands or unimproved grasslands, and 0.25 for forests), reflecting differences in residue quality and lignin content [7]. This procedure yielded coherent, literature-based estimates of annual carbon inputs (t C ha⁻¹ yr⁻¹), which were subsequently used to drive RothC simulations of soil organic carbon dynamics under different land-use scenarios.

Climatic information on precipitation, temperature, and monthly evaporation in situ was obtained from the Aznalcázar meteorological station, located at 37º09'06'' N, 06º16'24'' W, at an elevation of 2 metres above sea level. This database provided continuous records from October 2000 to February 2025. In parallel, to obtain climate records of precipitation, temperature, and evaporation derived from remote sensing, a script was implemented in the Google Earth Engine (GEE) processing environment. The ERA5-Land dataset [36] was employed, a high-resolution climate reanalysis product developed by the European Centre for Medium-Range Weather Forecasts (ECMWF) within the framework of the Copernicus Climate Change Service (C3S). This dataset is available from 1950 onwards, with a spatial resolution of approximately 9 km and a monthly temporal frequency. The selected variables were extracted and processed through GEE for the study period.

Finally, since radiocarbon measurements and SOC stock pools were not available, the RothC model was initialised through a spin-up procedure using the SoilR implementation. First, the monthly environmental modifier (ξ) was calculated from temperature (fT), water balance (fW; soil depth of 20 cm, pClay), and a constant cover factor (fC = 0.6), under the long-term mean climate. Subsequently, a function was used to derive an equilibrium allocation of the observed initial SOC stock across the DPM, RPM, BIO, HUM, and IOM compartments, assuming the following litter quality ratios for forest, shrubland, and grassland, respectively: DR = 0.25, DR = 0.67, and DR = 1.44. The IOM compartment was not prescribed using a pedotransfer function but emerged from the equilibrium solution. This procedure also provided a baseline carbon input (≈ 2.20 Mg C ha⁻¹ yr⁻¹), which was then used as a reference for transient simulations with time-varying inputs derived from the litter table.

2.4. Model Spin-Up Initialization

Before running the RothC simulations, a spin-up procedure was applied to initialise the distribution of SOC among the model pools. This step was necessary because direct measurements of the different RothC carbon compartments were not available for the study area. The spin-up was implemented using the SoilR package in R, following the RothC model structure. For each land-cover type, the observed initial SOC stock from the baseline sampling period was used as the target SOC value. The model was then run under long-term mean climatic conditions until an equilibrium state was reached. Monthly environmental modifiers were calculated from temperature, soil moisture and vegetation cover effects. The soil moisture modifier was estimated using monthly precipitation and evapotranspiration data, soil depth and clay content, while a constant vegetation cover factor was applied to represent the protective effect of plant cover on decomposition. The quality of organic inputs was defined through land-cover-specific DPM/RPM ratios. A ratio of 0.25 was used for forested areas, 0.67 for shrubland, and 1.44 for grassland, reflecting differences in residue decomposability and lignin content. These ratios determined the allocation of fresh organic matter inputs between the fast-cycling and more resistant plant material pools.

The spin-up procedure allowed the model to derive an internally consistent initial allocation of SOC across the RothC pools. This ensured that the baseline SOC stock used in the simulations was compatible with local climate, soil texture and land-cover conditions. The resulting pool distribution was subsequently used as the initial condition for the transient simulations covering the 2000–2024 period. In this way, the spin-up reduced the uncertainty associated with unknown initial pool sizes and provided a stable starting point for comparing the different boundary condition scenarios.

2.5. Sublandscape Analysis Through Boundary Conditions Scenarios (BC0-BC7)

To evaluate the performance of the RothC model in simulating soil organic carbon (SOC) dynamics in the sampling plots of Guadiamar Green Corridor, a boundary condition-based approach was implemented (Table 3) across the three main land-use types considered in this study: grassland, shrubland and riparian forest. Boundary conditions (BCs) were defined as alternative model-input configurations in which specific field-derived variables were replaced by values obtained from remote sensing products or global soil databases. This approach allowed us to assess the extent to which non-field-based data sources can be used to support SOC modelling and monitoring, reporting and verification (MRV) methodologies under different levels of data availability.

The analysis started with the reference scenario, BC0, in which all RothC input variables were derived from field measurements or local meteorological records, with the exception of carbon inputs. Since direct records of annual carbon inputs were not available for the whole study period, this variable was estimated from vegetation dynamics using NDVI-derived proxies, as described in the previous section. BC0 was therefore considered the closest approximation to an in situ-based modelling configuration and was used as the baseline scenario against which the remaining boundary conditions were compared. BC1 was designed to evaluate the effect of replacing the initial SOC stock with values derived from SoilGrids. In this scenario, the baseline SOC used to initialise the RothC model was obtained from global soil information rather than from field observations, while the remaining variables were kept consistent with the reference configuration. This scenario was included to assess whether globally available SOC datasets can provide a reliable alternative for RothC initialisation in areas where field measurements are scarce or unavailable. BC2 assessed the influence of replacing clay content with SoilsGrids values. Clay is a key parameter in RothC because it influences the stabilisation of organic carbon and the partitioning of decomposed material between CO₂, microbial biomass and humified organic matter. Therefore, this scenario allowed us to test the sensitivity of SOC dynamics to modelled soil texture information. BC3 to BC6 focused on the substitution of climatic variables. In BC3, precipitation data were replaced by remote sensing or reanalysis-derived estimates; in BC4, air temperature was substituted; and in BC5, evapotranspiration was replaced. BC6 combined the substitution of all climatic variables, providing an integrated assessment of the effect of using remote sensing-based climate inputs instead of local meteorological station records. These scenarios were used to evaluate the suitability of open-access climate datasets for driving RothC simulations in contexts where long-term local climatic records are incomplete or unavailable. Finally, BC7 evaluated the effect of replacing all the variables with the remote sensing derived variables.

2.6. Landscape Analysis

2.6.1. Data Acquisition and Spatiotemporal Structuring

Baseline soil organic carbon stocks $\left(\mathrm{S}\mathrm{O}{\mathrm{C}}_0\right)$ and soil clay content (%) were obtained from the SoilGrids and OpenLandMap global datasets at a 250 m spatial resolution. Monthly hydrometeorological forcing data (air temperature, precipitation, and potential evapotranspiration) were extracted from the ECMWF ERA5-Land reanalysis dataset for the 2000–2024 period. Annual soil carbon inputs $\left(C_{\mathrm{input}}\right)$ were estimated using Net Primary Productivity (NPP) products from the MODIS sensor (MOD17A3HGF). All spatial layers were reprojected to the UTM Zone 31N coordinate reference system (EPSG:32631) and bilinearly resampled to the 250 m reference grid. To optimize computational efficiency, the spatial domain was vectorized into a 2D spatiotemporal matrix $\left(N_{\mathrm{pixels}}\times 300\mathrm{months}\right)$. Because the satellite-derived NPP data were annual, carbon inputs were temporally disaggregated into monthly intervals using an empirical weighting vector derived from Mediterranean ecosystem litterfall dynamics $(w=\left[3.8,6.9,4.7,2.7,5.6,6.6,5.7,11.7,12.3,24.7,10.8,4.2\right]$ from January to December, capturing the characteristic autumn litterfall peak.

2.6.2. Model Parameterization, Spin-Up, and Transient Simulation

The RothC model was parameterized with a soil depth profile of 20 cm, a land cover factor $\left(f_c\right)$ of 0.6, zero external farmyard manure additions $\left(FYM=0\right)$ and a decomposable-to-resistant plant material (DPM/RPM) ratio of 0.2 (representing 17% decomposable and 83% resistant biomass). Default specific decomposition rates at 20 °C were retained (10.0, 0.3, 0.66, and 0.02 ${\mathrm{yr}}^{-1}$ for the DPM, RPM, BIO, and HUM pools, respectively). To establish thermodynamic equilibrium between the active carbon pools and local edaphoclimatic conditions, an analytical spin-up procedure was performed. Steady-state initial carbon compartments $(C_0$) were resolved via matrix inversion using the SoilR R package, driven by the 2000–2010 baseline climatology. Following spin-up equilibrium, pixel-by-pixel transient simulations were executed dynamically across the 2000–2024 historical period. At each monthly time step, carbon pool stocks were updated based on local temperature and soil moisture rate-modifying factors. From the resulting spatiotemporal matrix, annual SOC stock rasters were extracted for each modeled year, alongside a final spatial map of net cumulative SOC stock changes $\left(\Delta SOC=SO{C}_{2024}-SO{C}_{2000}\right)$, expressed in $\left(\mathrm{t}\mathrm{C}\cdot {\mathrm{ha}}^{-1}\right)$.

2.6.3. Spatial Uncertainty Propagation and Raster Outputs

To account for landscape heterogeneity and input data uncertainty, a spatially explicit bounding scenario framework was implemented. For each valid grid cell, upper and lower simulation bounds were calculated by simultaneously perturbing annual carbon inputs $\left(C_{\mathrm{input}}\right)$ and the decomposition rate modifier $\left(DR\right)$ by $20$ This yielded a conservative lower-bound scenario $\left(C_{\mathrm{input}}\times 0.8,DR\times 0.8\right)$ and an optimistic upper-bound scenario $\left(C_{\mathrm{input}}\times 1.2,DR\times 1.2\right).$ Based on these bounding runs, spatial uncertainty metrics were derived for the final year (2024), mapping the 5th and 95th carbon stock percentiles $\left(SO{C}_{\mathrm{p}05}andSO{C}_{\mathrm{p}95}\right)$, the absolute uncertainty range $\left(SO{C}_{\mathrm{p}95}-SO{C}_{\mathrm{p}05}\right)$ , and the relative spatial uncertainty via the Prediction Interval Ratio $\left(PIR=\left[SO{C}_{\mathrm{p}95}-SO{C}_{\mathrm{p}05}\right]/SO{C}_{2024}\right)$. The workflow yielded a high-resolution GeoTIFF database comprising monthly and annual SOC time-series raster stacks, baseline and final SOC inventories, the cumulative $\Delta SOC$ balance map, and the comprehensive suite of spatial uncertainty layers.

3. Results

3.1. RothC Model Performance and Evaluation of Sublandscape Results

The performance of the RothC model varied substantially among land-use types and boundary condition scenarios (Table 4 and Table 5). Overall, the results indicate that model accuracy was strongly influenced by both the vegetation cover and the type of input data used to initialise and drive the simulations.

Among the three land-use classes, woodland showed the best model performance. The lowest errors were obtained under BC6 and BC5, with RMSE values reaching approximately 7.26 Mg C ha⁻¹ in the best-performing scenario. This was also supported by a Nash–Sutcliffe efficiency coefficient of 0.65 and a Willmott agreement index of 0.88, indicating a robust agreement between observed and modelled SOC dynamics. These results suggest that RothC was able to reproduce SOC evolution more accurately in woodland areas, probably because carbon inputs and soil conditions were more consistent with the assumptions of the model. Grassland showed an intermediate model performance, with RMSE values generally around 14 Mg C ha⁻¹. Although the model was able to capture the general SOC trend, the lower efficiency values and agreement index indicate a higher dispersion between observed and simulated values than in woodland. This suggests that grassland SOC dynamics were moderately represented by RothC, but with greater uncertainty in the estimation of carbon inputs and temporal variability.

Shrubland presented the highest variability among boundary condition scenarios. While some configurations, particularly BC6, achieved relatively good performance, with an RMSE of 14.86 Mg C ha⁻¹, other scenarios produced much larger errors. In particular, BC1 generated RMSE values above 30 Mg C ha⁻¹, indicating a strong sensitivity of shrubland simulations to changes in boundary conditions. This behaviour suggests that shrubland SOC dynamics are more difficult to reproduce, probably due to greater spatial heterogeneity, variable vegetation recovery and uncertainty in the estimation of initial SOC and carbon inputs.

The comparison among boundary condition scenarios showed that BC6 was the best-performing configuration across all land-use types (Table 6). This scenario presented the lowest RMSE values and the highest model efficiency values, indicating the most consistent agreement between observed and simulated SOC stocks. BC5 also showed similar results, particularly at landscape scale, and can therefore be considered a solid alternative configuration for SOC modelling in the study area.

In contrast, BC1 showed the poorest performance, especially in woodland and shrubland. This scenario produced the largest discrepancies between observed and modelled SOC values and showed a marked negative PBIAS, indicating a systematic underestimation of SOC. This confirms that the model was particularly sensitive to the replacement of the initial SOC stock, since this parameter directly controls the allocation of carbon into the RothC pools during model initialisation.

The PBIAS analysis showed that BC0, BC2, BC3, BC4 and BC7 tended to slightly overestimate SOC stocks in woodland and shrubland, whereas BC5 and BC6 produced values closer to zero, indicating a more balanced representation of SOC dynamics. These results reinforce the suitability of BC6 for long-term projections in the Guadiamar Green Corridor.

Overall, these results indicate that the BC6 configuration provided the most reliable representation of the observed SOC dynamics in the Guadiamar Green Corridor. The model performed particularly well in woodland areas, while grassland and shrubland showed greater uncertainty. The poor performance of BC1 highlights the importance of accurately estimating initial SOC stocks before using remote sensing-derived inputs in process-based SOC models.

3.2. Boundary Conditions

Regarding the analysis of soil organic carbon (SOC) dynamics for each land use and boundary condition, we obtained the results shown in Figure 3. Overall, the simulated SOC trends over the 25-year period varied substantially depending on the land-use type and the boundary condition scenario applied. Among the three land-use classes evaluated, woodland showed the closest alignment between the simulated dynamics and the observed field data. For this land cover, the curves for scenarios such as BC6 and BC5 successfully mirrored the overall SOC trends, showing less deviation from the measured points throughout the simulation period. Grasslands demonstrated an intermediate performance. While the model captured the general temporal trend of SOC dynamics in this land use, the graph shows a slightly higher visual dispersion between the simulated lines and the observed data points compared to woodlands. Shrublands presented the highest variability among the boundary condition scenarios. While some configurations, particularly BC6, traced the observed SOC changes relatively well, other scenarios like BC1 exhibited large discrepancies, failing to capture the actual dynamics accurately.

When comparing the different boundary condition scenarios across all land uses, the curves for BC6 consistently showed the best visual fit to the observed SOC changes. Scenario BC5 also displayed very similar and consistent trends. Conversely, the curve for scenario BC1 showed the poorest performance; the graph illustrates that this scenario systematically underestimated the SOC stocks from the beginning of the simulation, presenting the largest gap between the modelled dynamics and the actual field measurements across the different land uses.

3.3. Landscape Results

Across the 1,008 valid 250 m grid cells (Figure 4), the mean soil organic carbon (SOC) stock in the Guadiamar basin increased from 22.68 t C ha⁻¹ in 2000 to a projected 25.45 t C ha⁻¹ in 2024. This represents a net regional carbon gain $\left(\Delta SOC\right)$ of $2.77\pm 2.62$ t C ha⁻¹. Approximately 85.02% of the landscape (857 pixels) accumulated carbon, with a mean gain of $+3.59$ t C ha⁻¹ within these areas and peak sequestration values reaching up to $+10.21$ t C ha⁻¹. Soil carbon losses occurred across 14.98% of the landscape (151 pixels), averaging $-1.86$ t C ha⁻¹ with maximum localized losses reaching $-6.42$ t C ha⁻¹.

The observed SOC accumulation was strongly and positively correlated with annual MODIS-derived organic inputs $\left(r=0.743\right)$ . Specifically, net carbon sink areas received a mean input of 3.48 Mg C ha⁻¹ yr⁻¹, compared to 2.89 Mg C ha⁻¹ yr⁻¹ in carbon source areas.Furthermore, a pronounced soil carbon saturation effect was identified; $\Delta SOC$ exhibited a strong negative correlation with the initial baseline stock from the year 2000 $\left(r=-0.797\right)$ Grid cells that acted as carbon sources started with significantly higher baseline stocks (26.34 t C ha⁻¹) than those functioning as sinks (22.03 t C ha⁻¹).

Spatially explicit bounding scenario analysis $(\pm 20$ perturbation in both carbon inputs and decomposition rate modifiers) yielded a mean landscape-wide relative uncertainty of 11.28% as showed in Figure 5 (ranging narrowly between 7% and 14%). The robust nature of the simulation is highlighted by the fact that even under the conservative 5th percentile scenario $\left(SO{C}_{\mathrm{p}05}\right)$, the net regional carbon balance remained distinctly positive. Propagation of input parameter uncertainty through the bounding scenario framework $(\pm 20$ perturbations in annual carbon inputs and decay rate modifiers) demonstrated high spatial robustness across the regional modeling domain.

Across all 1,008 valid grid cells, the untransformed PIR index yielded a regional mean ratio of 0.226 (or 22.56% when expressed as percentage of the central inventory). Spatially, the PIR ratio remained remarkably stable across the basin, with the 50% interquartile range constrained between 0.213 (21.29%) and 0.238 (23.77%), and absolute extremes ranging from a minimum of 0.156 (15.55%) in highly stable edaphoclimatic sectors to a maximum of 0.280 (28.03%) in zones highly sensitive to input perturbations. Notably, spatial PIR values exhibited a tight positive correlation with $\Delta SOC$ magnitude $\left(r=+0.991\right)$ , confirming that landscape sectors undergoing the most vigorous carbon accumulation or rapid depletion naturally carry proportionally higher absolute sensitivity to boundary parameter perturbations.

4. Discussion

4.1. Sub-Landscape

Among the evaluated land-cover types, riparian forests exhibited the highest modelling consistency and predictive accuracy. Under the best-performing scenario (BC6), woodland simulations achieved an exceptional statistical alignment with observed field inventories . This robust performance can be attributed to the ecophysiological stability of Mediterranean riparian forests in the Guadiamar Green Corridor. In established woodland formations, deep rooting systems, dense canopy cover, and a persistent litter layer create a highly buffered microclimate that insulates the topsoil from extreme temperature spikes and rapid desiccation cycles [34]. Furthermore, the parameterisation of organic matter partitioning in RothC accurately reflected the high proportion of recalcitrant, lignin-rich structural tissues characteristic of woody litterfall. Consequently, the decomposition kinetics simulated by RothC under standard rate constants closely mirrored the natural biogeochemical equilibrium of these restored afforested plots.

In contrast, grassland and shrubland formations presented intermediate to high modelling uncertainty, accompanied by lower efficiency coefficients). In Mediterranean grasslands, carbon turnover is inherently dynamic and subject to pronounced interannual climatic fluctuations. The high allocation of net primary productivity to fine roots and easily decomposable herbaceous tissues makes grassland topsoils highly responsive to seasonal precipitation pulses and summer drought dormancy [32]. However, the interpretation of grassland model performance requires additional caution because this land-cover class was represented by a single long-term monitoring plot. Therefore, the grassland results provide useful information on model behaviour under herbaceous vegetation within the Guadiamar restoration context, but they should not be interpreted as a fully replicated assessment of grassland SOC dynamics at landscape scale. In shrublands, the observed temporal dispersion stems from the extreme spatial heterogeneity of vegetation colonization and structural patchiness following the post-mining restoration. Patchy shrub canopies generate highly localized islands of fertility separated by bare or sparsely vegetated soil, resulting in high micro-scale variability in field SOC measurements that a 1D point-based simulation cannot fully capture without high-resolution spatial aggregation.

The progressive substitution of in situ observations with remote sensing and reanalysis datasets (the boundary condition framework, BC0–BC7) provided critical insights into the operational feasibility and sensitivity of process-based soil organic carbon (SOC) modelling in restoration contexts. A standout finding of this study is the remarkable performance of scenario BC6, where local meteorological station records of temperature, precipitation, and evapotranspiration were entirely replaced by the ERA5-Land reanalysis dataset. Across all evaluated land uses—most notably in woodland formations $\left(\mathrm{RMSE}=\mathbf{7.26}\mathrm{Mg}\mathrm{C}\cdot {\mathrm{ha}}^{-1},\mathrm{EF}=\mathbf{0.65}\right)$—the BC6 configuration matched or even slightly exceeded the predictive accuracy of the in-situ reference scenario (BC0).

A critical biogeochemical anomaly identified in the sublandscape simulation is the widespread, systematically decreasing trend in soil organic carbon (SOC) stocks observed across several monitoring plots over the 25-year timeframe. Rather than indicating an ecosystem degradation process, this downward trajectory represents a transitional relaxation phase toward a natural thermodynamic equilibrium. The initial soil inventories measured in 2000 were artificially inflated due to the massive application of organic-clay amendments implemented as an emergency remediation strategy to immobilize heavy metals and reconstruct the topsoil profile immediately after the 1998 Aznalcóllar mine spill [18]. These baseline inventories far exceed the natural carrying capacity of typical Mediterranean Fluvisols under dryland conditions, which baseline models and empirical regional studies estimate to range between 20 and 35 [39,40]. Once the active external organic inputs ceased, the highly active microbial pools simulated in RothC rapidly mineralized the labile fractions of these exogenous amendments. Crucially, only the riparian forest (woodland) formations possessed high net primary productivity (NPP), stable litterfall inputs, and buffered microclimatic moisture conditions necessary to physically protect and sustain these elevated carbon stocks . In contrast, the water-limited grasslands and shrublands, subjected to intense Mediterranean summer drought and lower organic return rates, could not support the amended carbon inventories, triggering a steep transient decay toward their true biogeochemical carrying capacity.

This initial anthropogenic distortion was further compounded by extreme micro-spatial heterogeneity within the experimental plots, as evidenced by the exceptionally high standard deviations (STD) recorded among sampling replicates during the baseline campaign in the year 2000. In these early post-remediation stages, the spatial variance reached critical magnitudes relative to the plot means. Most notably in woodland plots like, where the standard deviation represented nearly 40% to 45% of the total carbon stock. This massive spatial noise serves as a direct diagnostic signature of the highly uneven, mechanically patchy application of composted amendments and red clays during the emergency soil reconstruction phase [17]. Over the 25-year chronosequence, however, this early spatial heterogeneity experienced a highly pronounced natural stabilization and homogenization. By 2014 and 2024, the standard deviations between replicates dropped dramatically across almost all monitoring plots. This long-term spatial stabilization reflects the gradual biogeochemical integration of the exogenously added organic matter through continuous root expansion, macro-invertebrate bioturbation, and localized microbial mineralization of carbon hotspots [20]. These findings confirm that while empirical baseline inventories in reconstructed soils are highly vulnerable to micro-scale spatial noise, natural soil development processes act as spatial homogenizers over multi-decadal scales, gradually reconciling local field inventories with regional baseline projections

These contrasting outcomes underscore a key methodological principle for large-scale MRV frameworks: remote sensing products and global databases cannot be applied as universal substitutes for all model inputs. While continuous meteorological drivers and relatively stable soil physical parameters, such as clay content, may be reliably replaced by open-access reanalysis or gridded datasets, the initial SOC baseline (${\mathit{S}\mathit{O}\mathit{C}}_{\mathbf{0}}$) remains a critical parameter that requires empirical ground-truthing. This is particularly important in restored or heavily disturbed systems, where historical management legacies, organic amendments, soil reconstruction practices, and local contamination-remediation actions may generate baseline SOC stocks that strongly deviate from the natural biogeochemical carrying capacity of the system. In such cases, process-based models such as RothC may reproduce the general direction of SOC change, but they cannot fully recreate SOC dynamics under all situations, especially when the initial carbon stock is far above what the soil, vegetation, and bioclimatic conditions can sustain over the long term. Therefore, accurate baseline sampling is essential not only for model initialisation, but also for interpreting whether subsequent SOC losses reflect model failure, ecosystem degradation, or a natural adjustment towards a new equilibrium state. These findings highlight the need to combine field observations, remote sensing products, and process-based simulations within a carefully validated framework, rather than treating spatial datasets as direct replacements for local soil information.

4.2. Landscape

The landscape-scale simulation provides a complementary perspective to the plot-based boundary condition analysis. However, it is important to emphasise that this approach was not based on SOC measurements. Instead, all input variables were derived from remote sensing products, climate reanalysis data, and global soil databases. In particular, the SOC values used to initialise the year 2000 baseline were obtained from SoilGrids, representing an attempt to reconstruct the spatial baseline condition of the Guadiamar Green Corridor using globally available soil information. Therefore, the landscape analysis should be interpreted as a spatially explicit modelling exercise based on open-access datasets, rather than as a direct reconstruction of measured SOC stocks across the restoration area. Under this configuration, the model simulated an overall increase in SOC stocks between 2000 and 2024, with most of the landscape acting as a carbon sink. This result contrasts with some of the plot-scale observations, where several sampling plots showed a decline in SOC over the same period. This apparent discrepancy is mainly explained by the different nature of the baseline used in each analysis. At plot scale, the 2000 baseline reflected field-measured SOC stocks shortly after remediation, when large amounts of organic and clay-rich amendments had been applied to the most contaminated areas. These initial SOC values were therefore strongly influenced by local restoration practices and, in several plots, were considerably higher than what the soil–vegetation–climate system could sustain over the long term. By contrast, the SoilGrids-derived baseline used in the landscape simulation represents a modelled regional estimate of SOC, which is more likely to reflect broader edaphoclimatic conditions than the highly localised effects of post-spill amendment applications.

This distinction is critical for interpreting the results. The positive SOC balance obtained at landscape scale does not necessarily contradict the SOC decline observed in some field plots. Rather, it suggests that, when initial SOC stocks are defined from a global soil product and driven by remotely sensed productivity and reanalysis climate data, the Guadiamar landscape has the potential to function as a moderate carbon sink under current vegetation and climatic conditions. In other words, the landscape simulation represents the expected SOC trajectory under a harmonised, gridded baseline, whereas the sub-landscape analysis captures the legacy of site-specific remediation actions and the subsequent relaxation of artificially elevated SOC stocks towards a new equilibrium.

The use of SoilGrids as the baseline dataset also introduces an important methodological limitation. Global digital soil mapping products are designed to provide spatially continuous estimates over large areas, but they may not capture abrupt, localised or anthropogenic soil modifications, such as those associated with mining-spill remediation, topsoil removal, clay amendment application or organic amendment incorporation. This reinforces the conclusion obtained from the boundary condition analysis: global soil datasets can be useful for regional modelling where field data are unavailable, but their outputs must be interpreted cautiously in disturbed or reconstructed landscapes.

The strong relationship between simulated SOC accumulation and MODIS-derived carbon inputs further indicates that vegetation productivity was a major driver of landscape-scale SOC dynamics. MODIS GPP/NPP products are specifically designed to provide spatially continuous estimates of terrestrial vegetation productivity and are widely used for ecosystem status assessment, carbon-cycle analysis and environmental monitoring [42,43,44]. Areas receiving higher annual organic inputs tended to behave as carbon sinks, whereas areas with lower productivity or higher initial SOC stocks were more likely to show limited accumulation or carbon losses. This pattern is consistent with recent discussions on soil carbon saturation, which emphasise that SOC accumulation is constrained by the capacity of soils to stabilise additional carbon, and that the initial SOC status strongly influences the potential for further sequestration [45,46]. In this context, soils with lower initial SOC or a larger saturation deficit generally have greater potential to accumulate additional carbon, while soils already close to their bioclimatic or mineral-associated stabilisation capacity tend to accumulate carbon more slowly or may even lose carbon over time.

The uncertainty analysis showed that the simulated carbon sink remained relatively robust under perturbations of carbon inputs and decomposition modifiers. Nevertheless, this robustness should be understood within the assumptions of the gridded modelling framework. The relatively stable prediction interval ratio indicates internal consistency of the spatial RothC implementation, but it does not eliminate the uncertainty associated with the choice of baseline dataset. SoilGrids 2.0 provides globally consistent soil-property predictions at 250 m resolution with spatially explicit uncertainty, using machine-learning models trained on global soil observations and environmental covariates [47]. However, recent studies also highlight that global soil raster products may provide useful regional-scale information while still showing important local inaccuracies when evaluated against external soil datasets, especially in heterogeneous landscapes or areas affected by specific local processes [48]. In this sense, the largest source of uncertainty is not only the propagation of model parameters, but also the representativeness of the initial SOC map used to initialise the simulation. This limitation is particularly relevant in the Guadiamar Green Corridor, where the year 2000 baseline derived from SoilGrids represents a modelled regional approximation rather than a field-measured post-remediation SOC inventory.

Overall, the landscape analysis demonstrates the potential of combining SoilGrids, OpenLandMap, ERA5-Land and MODIS-derived productivity data to implement spatially explicit RothC simulations in areas where complete field datasets are unavailable. However, the results should be interpreted as a regional modelling scenario based on remote sensing and global baseline products, not as a direct substitute for field-based SOC monitoring. For restored and anthropogenically disturbed systems such as the Guadiamar Green Corridor, landscape-scale simulations are most valuable when used together with plot-scale observations, allowing the identification of discrepancies between modelled background SOC trajectories and local soil dynamics shaped by remediation history.

5. Conclusions

This study demonstrates the potential and limitations of using RothC to simulate long-term soil organic carbon dynamics in a Mediterranean restoration landscape affected by a major mining spill. The Guadiamar Green Corridor represents a particularly complex case study because its initial SOC conditions were strongly shaped by post-disaster remediation actions, including the application of organic and clay-rich amendments. As a result, the SOC baseline measured shortly after restoration did not necessarily represent a stable edaphoclimatic equilibrium, but rather a transient and highly heterogeneous post-remediation condition.

At the sub-landscape scale, RothC was able to reproduce the general direction of SOC change observed over the 25-year period, although its performance varied substantially among land-cover types and boundary condition scenarios. Riparian woodland showed the most consistent model performance, suggesting that its higher vegetation productivity, litter inputs and more buffered microclimatic conditions favoured a more predictable SOC trajectory. In contrast, grassland and shrubland areas showed greater uncertainty, reflecting stronger temporal variability, spatial heterogeneity and a higher sensitivity to the estimation of carbon inputs and initial SOC conditions.

The boundary condition analysis confirmed that not all RothC input variables have the same influence on model performance. The replacement of local meteorological records with ERA5-Land climate variables produced robust results, particularly under the BC6 configuration, indicating that reanalysis datasets can be a reliable alternative for driving process-based SOC simulations where long-term local climate records are unavailable. However, the substitution of the initial SOC baseline produced the largest deviations, especially in scenarios using globally modelled SOC estimates. This highlights that the initial SOC stock is not merely an input parameter, but a structural condition that controls pool allocation, decomposition trajectories and the subsequent interpretation of SOC change.

The results therefore show that remote sensing products and global environmental datasets can support SOC modelling, but they cannot be applied as universal replacements for field data. This is especially relevant in restored or anthropogenically disturbed systems, where local management history, soil reconstruction and amendment application may generate SOC stocks that differ markedly from the values expected under natural soil, vegetation and bioclimatic conditions. Under these circumstances, process-based models may capture broad tendencies, but they cannot be expected to fully reproduce all SOC dynamics when the baseline state is highly artificial or far from equilibrium.

At the landscape scale, the spatial RothC implementation provided a complementary regional perspective based entirely on gridded datasets. In this case, the year 2000 baseline was not derived from field measurements, but from SoilGrids-derived SOC estimates, while clay content, climate forcing and carbon inputs were obtained from OpenLandMap, ERA5-Land and MODIS-derived productivity data. The simulated regional SOC increase should therefore be interpreted as a remote sensing-based modelling scenario under a globally estimated baseline, rather than as a direct reconstruction of field-observed SOC changes across the Guadiamar Green Corridor.

The apparent contrast between the plot-scale decline in SOC and the landscape-scale carbon gain is mainly explained by the different baseline definitions. The field-based analysis captured the legacy of post-remediation amendments and the subsequent relaxation of artificially elevated SOC stocks towards a new equilibrium. In contrast, the landscape analysis represented a broader gridded baseline that was more closely aligned with regional edaphoclimatic conditions and therefore suggested a moderate carbon sink potential under current vegetation productivity. This distinction reinforces the importance of baseline selection as one of the most influential decisions in spatial SOC modelling.

The uncertainty analysis indicated that the spatial model was internally robust under perturbations of carbon inputs and decomposition modifiers. Nevertheless, the main source of uncertainty remains the representativeness of the initial SOC map used to initialise the simulations. For this reason, global soil products such as SoilGrids can be valuable for regional screening and modelling in data-scarce contexts, but their outputs should be validated against field observations whenever the study area has experienced strong local disturbance or restoration interventions.

Overall, this study supports the use of RothC, remote sensing and open-access environmental datasets as complementary tools for scalable SOC monitoring. However, it also demonstrates that robust SOC modelling in restored landscapes requires a hybrid framework in which field measurements are used to constrain and validate model baselines, while remote sensing and gridded datasets provide spatial continuity and temporal consistency. In Mediterranean restoration contexts such as the Guadiamar Green Corridor, this combined approach is essential to distinguish between real ecosystem degradation, natural SOC stabilisation after remediation, and model uncertainty arising from baseline assumptions.