Machine learning in space and time for modelling soil organic carbon change

Spatially resolved estimates of change in soil organic carbon (SOC) stocks are necessary for supporting national and international policies aimed at achieving land degradation neutrality and climate change mitigation. In this work we report on the development, implementation and application of a data‐driven, statistical method for mapping SOC stocks in space and time, using Argentina as a pilot. We used quantile regression forest machine learning to predict annual SOC stock at 0–30 cm depth at 250 m resolution for Argentina between 1982 and 2017. The model was calibrated using over 5,000 SOC stock values from the 36‐year time period and 35 environmental covariates. We preprocessed normalized difference vegetation index (NDVI) dynamic covariates using a temporal low‐pass filter to allow the SOC stock for a given year to depend on the NDVI of the current as well as preceding years. Predictions had modest temporal variation, with an average decrease for the entire country from 2.55 to 2.48 kg C m−2 over the 36‐year period (equivalent to a decline of 211 Gg C, 3.0% of the total 0–30 cm SOC stock in Argentina). The Pampa region had a larger estimated SOC stock decrease from 4.62 to 4.34 kg C m−2 (5.9%) during the same period. For the 2001–2015 period, predicted temporal variation was seven‐fold larger than that obtained using the Tier 1 approach of the Intergovernmental Panel on Climate Change and United Nations Convention to Combat Desertification. Prediction uncertainties turned out to be substantial, mainly due to the limited number and poor spatial and temporal distribution of the calibration data, and the limited explanatory power of the covariates. Cross‐validation confirmed that SOC stock prediction accuracy was limited, with a mean error of 0.03 kg C m−2 and a root mean squared error of 2.04 kg C m−2. In spite of the large uncertainties, this work showed that machine learning methods can be used for space–time SOC mapping and may yield valuable information to land managers and policymakers, provided that SOC observation density in space and time is sufficiently large.

temporal distribution of the calibration data, and the limited explanatory power of the covariates. Cross-validation confirmed that SOC stock prediction accuracy was limited, with a mean error of 0.03 kg C m −2 and a root mean squared error of 2.04 kg C m −2 . In spite of the large uncertainties, this work showed that machine learning methods can be used for space-time SOC mapping and may yield valuable information to land managers and policymakers, provided that SOC observation density in space and time is sufficiently large.

| INTRODUCTION
The monitoring, modelling and mapping of soil organic carbon (SOC) is important for many reasons. SOC, a proxy for soil organic matter content, is a major determinant of soil quality and soil fertility, SOC sequestration can contribute substantially to climate change mitigation (Batjes, 2019;Minasny et al., 2017), and SOC stock is one of three Land Degradation Neutrality indicators used by the United Nations Convention to Combat Desertification (UNCCD). Parties to the UNCCD have agreed to report on SOC stock trends at regular time intervals in Nationally Determined Contributions (UNCCD, 2019).
By building on recent advances in digital soil mapping (DSM), many SOC maps on national, continental and global scales have been produced during the past decade (e.g., Adhikari et al., 2014;FAO & ITPS, 2018;Hengl et al., 2017;Kempen et al., 2019;Mulder, Lacoste, Richerde-Forges, & Arrouays, 2016;Stoorvogel, Bakkenes, Temme, Batjes, & Ten Brink, 2017). However, the majority of these maps are static, whereas SOC is dynamic and SOC dynamics are of particular interest to carbon sequestration and land degradation studies. Thus, there is a clear need to extend spatial SOC mapping to space-time SOC mapping.
There are not that many studies on modelling and mapping SOC variation in space and time at the national to global scale, but several approaches can be distinguished. Perhaps the simplest approach is that used in the UNCCD-modified intergovernmental panel on climate change (IPCC) Tier 1 and Tier 2 methods (IPCC, 2006;Mattina et al., 2018). Here, one starts with a static, baseline SOC stock map for a reference year, and models the SOC change from the reference year onward by modification of the baseline map through multiplication with land-use change, management and input factors. This approach is commonly adopted by countries when reporting on carbon trends to the UNCCD. Although it is relatively simple, the disadvantages are that it is not obvious how to derive a baseline SOC map for a reference year, that management and input factors are not easily obtained at the country scale (and hence in practice often ignored), and that linking SOC change directly to land-use change is challenging. Land-use change in itself is difficult to model (Verburg et al., 2019), strongly influenced by the land-use classification system used, and there is considerable SOC variation within land-use classes (Lamichhane, Kumar, & Wilson, 2019).
Another approach to modelling and mapping SOC in space and time is to first derive a static, empirical DSM model in the usual way and model the SOC temporal variation by recognizing that some of the covariates of the model are in fact dynamic. SOC temporal variation may then be derived by submitting these dynamic covariates to the DSM model. Using this approach, Stockmann et al. (2015) derived a global topsoil SOC model that had land-cover as a main driver of SOC change. By running this model with moderate resolution imaging spectroradiometer (MODIS) land-cover maps for 2001 and 2009, they derived a global SOC change map between 2001 and 2009 at 1-km resolution. Gray and Bishop (2016) mapped SOC changes for New South Wales until 2070 by calibrating a decision tree DSM model using legacy soil data and climate, terrain, parent material, age and biota covariates, and running the model on different climate change projections. Yigini and Panagos (2016) used a similar space-for-time-substitution approach to map SOC stocks in Europe under future climate and land-cover changes. Huang, Hartemink, and Zhang (2019) mapped SOC stock for several years in the 1850-2002 period for the state of Wisconsin, USA. They first calibrated a random forest DSM model using SOC data from 1980 onwards and vegetation, climate and terrain covariates. By utilizing the fact that climate and land-use were important drivers of their model, substituting historic land-use and climate maps provided insight into SOC stock change over the past 150 years. Although the space-fortime-substitution approach is attractive because it allows derivation of SOC maps beyond the time period within which SOC observations are made, it builds on the assumption that the relationship between SOC and covariates can be extrapolated over time. Szatmári et al. (2019) mapped the SOC stock change between 1992 and 2010 for Hungary by fitting separate DSM models for both years. They benefited from the fact that in Hungary sampling sites were revisited and SOC observations at the sampling sites were available for both years. The advantage of this method is that no time extrapolations are made, but it only works if SOC data are available for both years and cannot predict SOC in years for which no or insufficient SOC data are available.
Approaches to capture SOC dynamics that do more justice to the underlying physical, chemical and biological processes make use of semi-mechanistic models, such as RothC and CENTURY (e.g., Abramoff et al., 2018;Gottschalk et al., 2012;Karunaratne, Bishop, Lessels, Baldock, & Odeh, 2015;Lugato, Panagos, Bampa, Jones, & Montanarella, 2014;Smith et al., 2005;Woolf & Lehmann, 2019). Alternatively, Earth System Models have also been used (e.g., Todd-Brown et al., 2014). The semi-mechanistic approach is particularly attractive when the goal of modelling is not just prediction but also understanding, because it can integrate existing knowledge on soil processes and account for different soil carbon pools, but application on national to global scales is challenging. Different processes are dominant in different parts of the world and must all be represented in a global model. Model parametrization is difficult because it is not realistic to assume that parameters are spatially invariant in case of large-scale applications. Moreover, boundary conditions and driving forces of semimechanistic models are often poorly known. Computational requirements may also be limiting for large-scale applications, although this is not a fundamental problem given the rapid developments in geo-computation. This paper has three main objectives. First, we extend on existing space-time SOC mapping approaches by developing a fully data-driven, machine learning model that is calibrated with SOC observations that cover the entire time period for which predictions are to be made and that uses static and dynamic covariates from that same period. We decrease the risk of unwarranted extrapolations by requiring that the same covariates are used for calibration and prediction. We also require that the covariates are publicly available for the globe and have sufficiently fine spatial and temporal resolution, thus allowing future applications of this work to map any region on the globe where temporal soil data are available. Second, we test the performance of the model for a pilot area and time period. We use the space-time topsoil SOC stock in Argentina between 1982 and 2017 as a case study. We analyse results by evaluating spatial patterns and temporal trends, the latter both for the country as a whole as well as for regions within Argentina. We also quantify the accuracy of the model predictions, both using 10-fold cross-validation and by deriving 90% prediction interval maps for all years. Third, we compare the results of our model for the pilot area with those obtained with the UNCCD-modified IPCC Tier 1 approach and interpret differences.

| Space-time digital soil mapping methodology
The modelling approach used is machine learning for DSM . Specifically, we used the quantile regression forest (QRF) algorithm (Meinshausen, 2006;Szatmári et al., 2019;Vaysse & Lagacherie, 2017) to establish a relation between topsoil SOC stock and environmental covariates. The main difference with purely spatial DSM is that the regression matrix used for model calibration is derived from a space-time instead of a spatial overlay of SOC observations on covariates. This requires that the year of soil sampling is known for all SOC observations and that dynamic covariates are available for all years in which SOC observations are made. Once the regression matrix is derived, model calibration, cross-validation and prediction are carried out in the usual way (e.g., Hengl et al., 2017). Specific details of these steps, such as model selection, choice of hyperparameters and computational aspects, are provided after the soil data and covariates of the pilot study have been introduced.

| Argentina terrain, climate, land-use and soils
Argentina has a land surface area of 2.78 million km 2 (about one-sixth of South America). Although the country is known for its spacious plains with fertile soils and humid climate, these conditions represent less than onethird of its surface area. The remaining space is dominated by arid and semi-arid conditions. The climate is predominantly temperate, although it is subtropical in the north and subpolar in the far south. Argentina's orography is characterized by the presence of mountains in the west and plains in the east, with altitudes decreasing generally from west to east (Gardi et al., 2014). The main land-cover types are natural and semi-natural terrestrial vegetation (60%), cultivated and managed lands (20%), and natural or semi-natural aquatic or regularly flooded vegetation (8%). According to the Argentinian Ministry of Environment and Sustainable Development, more than 30,000 km 2 of forest have been lost between 2007 and 2016. Argentina is a country of marked geological, geomorphological and climatic contrasts, and thus presents a wide variety of soils. Three USDA Soil Taxonomy soil orders (Mollisols, Entisols and Aridisols) cover almost three-quarters of the country (Supporting Information, Figure S1). Less extensive soil orders are Alfisols, Inceptisols, Andisols, Vertisols, Histosols and Ultisols.

| SOC concentration data for Argentina
The soil dataset compiled for this project contains 18,768 observations on SOC concentration from 5,480 soil profiles throughout Argentina. The data were collected from 1955 onwards (Supporting Information, Figure S2). The highest observation density is between 1970 and 1990, corresponding with the national soil survey plan of the country. Between 2005 and 2015, most data were collected by soil fertility projects. The majority of observations (about 68%) are from the Argentinian SISINTA soil database (Olmedo, Rodriguez, & Angelini, 2017), but to increase sample size additional smaller datasets were added. All observations were subjected to quality and consistency checks and conversion factors applied to harmonize SOC concentration data to the Walkley-Black method (Richter, Massen, & Mizuno, 1973).
Most of the datasets that were merged for this study correspond with soil observations down the profile at multiple depths, sometimes using fixed intervals, sometimes by soil horizon. For all observations, the top and bottom of the sampling depth interval were recorded. Generally, there are more SOC observations for the upper soil layers than for the deeper soil layers. Conversion to 0-30 cm, 30-60 cm and >60 cm depth intervals was achieved by weighted averaging of SOC observations per depth interval, assuming constant SOC concentration within each sampling interval. Histograms of the SOC concentration for the three depth intervals are provided in the Supporting Information ( Figure S3). Note that in this study, we only used the topsoil SOC (0-30 cm).

| Conversion to SOC stock
Topsoil SOC stock values at sampling sites were derived from observations of SOC concentration, bulk density and coarse fragments as follows: where BD is bulk density and CRF the proportion of coarse fragments. Here, all variables refer to the top 0-30 cm of the soil. The proportion of coarse fragments at each sampling location was derived from the betaSoilGrids2019 product (SoilGrids, 2019). Bulk density was measured for only 940 out of the total of 18,768 soil samples. We "filled the gaps" using a random forest pedotransfer function (PTF) that predicts BD from other soil properties and environmental covariates (Ramcharan, Hengl, Beaudette, & Wills, 2017). We calibrated the PTF using a randomly selected subset ( Figure 1 also reports the mean error (ME) and the concordance correlation coefficient (CCC) between predictions and independent observations. Note that the PTF model produced biased predictions in Argentina, which is not unexpected because we applied a PTF model calibrated on a global dataset to a local validation set.
Using the PTF, we expanded the bulk density dataset to all 18,768 soil samples of the Argentina soil dataset. Thus, only 940 of these are real bulk density observations, whereas all other bulk density values are PTF predictions. Although the predictions are only proxies of the "true" bulk density, we ignored PTF errors and treated all 18,768 BD values as if "error free". After conversion to BD values for the 0-30 cm depth interval using weighted averaging, 0-30 cm SOC stock observations were derived using Equation 1. Supporting Information Figure S4 shows a histogram of the topsoil SOC stock data. The total number of 0-30 cm SOC stock values available for calibration of the SOC stock machine learning model was 5,114.

| Covariates for Argentina
We only used publicly available global covariates to facilitate future extension of the model from the pilot area to the globe. The main source of static covariates was the database used by SoilGrids  and included a digital elevation model and derivatives, land-cover, long-term average climate variables and lithology. In particular we used MODIS products from which a variety of vegetation indices were derived from spectral bands, together with information about primary productivity and evapotranspiration.
Dynamic vegetation indices (i.e., normalized difference vegetation index (NDVI)) were derived from advanced veryhigh resolution radiometer (AVHRR) images, globally available from 1982 onwards. These indices were available as daily values and were aggregated to seasonal values for each year. These seasonal covariates were further processed by generating weighted averages over multiple years, going back in time. We did this so that the SOC stock for a given year need not only depend on the NDVI for the seasons in that same year, but might also depend on seasonal NDVI values in previous years, because vegetation change has a delayed effect on SOC change. For this we used an exponential decay function (see Figure 2). All NDVI covariates associated with three values of the exponential decay parameter a were included in the machine learning algorithm (i.e., a = 0.6, a = 0.8 and a = 0.9). Application of the exponential decay function requires that covariates are needed for years before the starting year of the prediction time series. In order to not shorten the time series of the SOC prediction maps, we assumed that the dynamic covariates for the years before the starting year were equal to the covariates of the starting year. To reduce computing time, we also applied a threshold of 0.1, meaning that we only included contributions from the past for which the exponential decay function is above the threshold.
The original covariates varied in spatial resolution from 100 m to 10 km and were brought to a common resolution of 250 m and re-projected to an equal-area projection to avoid areal distortion (de Sousa, Poggio, & Kempen, 2019). They were also brought to a common extent by clipping to the pilot area (i.e., a gridded map of Argentina from which all non-soil areas were filtered out).

| Model selection, cross-validation and prediction
Recursive feature elimination analysis (RFE, Guyon, Weston, Barnhill, & Vapnik, 2002) was used to select the best performing subset of static covariates. The RFE procedure is similar to backward regression, in that it starts with the maximum number of covariates and iteratively removes the weakest explanatory variable until a F I G U R E 1 Density scatter plots of predicted against observed bulk density for a validation subset of the WoSIS database (left) and independent Argentinian validation data (right). Solid line is the 1:1 line. ME, mean error; RMSE, root mean squared error; CCC, concordance correlation coefficient [Color figure can be viewed at wileyonlinelibrary.com] specified number of covariates is reached. By doing this for all possible numbers of covariates and plotting performance indices computed on "out-of-bag" observations, the optimal number of covariates to be included in modelling is derived. Next, we extended the set of covariates by adding all 12 dynamic NDVI covariates (i.e., all combinations of four seasons and three exponential decay parameters). The QRF hyperparameters mtry and ntree were set to default values, that is, the square root of the number of covariates and 500, respectively.
The performance of the final model was evaluated using 10-fold cross-validation. Thus, the dataset was randomly split into ten equally sized folds, and nine of these were used to calibrate the QRF model and predict the SOC stock for the remaining fold. This procedure was carried out 10 times, each time setting aside a different fold. We plotted the predictions against the independent observations and computed the ME, RMSE and CCC.
Predictions were made for the whole of Argentina at 250 m resolution for each year in the 1982-2017 period. Note that because we used QRF, predictions could be made for different quantiles. We used the 0.5-quantile (i.e., the median) as a prediction of the SOC stock and the 0.05-and 0.95-quantiles as the lower and upper limits of a 90% prediction interval, respectively. We also spatially aggregated the maps of the median SOC stock for the entire country and for specific regions within the country, and plotted the spatial averages as time series to analyse temporal SOC stock trends. Note that we could not derive the uncertainty (i.e., the 0.05-and 0.95-quantiles) associated with these spatial aggregates because this requires the spatial autocorrelation of prediction errors (Kros et al., 2012;Plaza et al., 2018), which is not computed by QRF.

| Implementation and computational aspects
We used GRASS GIS (GRASS, 2019) for covariates preparation, data storage and tiling of predictions. The modelling was carried out with R software (R, 2019), in particular the ranger package (Wright & Ziegler, 2017). The procedure was implemented in parallel using a High Performance Computing facility of Wageningen University. The predictions were parallelized using a tiling approach. The full workflow (model fitting, cross-validation and prediction on a 250 m grid for each of the 36 years) took approximately 3,000 CPU hr.

| Comparison with UNCCD-modified IPCC Tier 1 results
Maps of the 0-30 cm SOC stock change for Argentina between 2001 and 2015 were also derived separately using the UNCCD-modified IPCC approach. This approach was reviewed in the Introduction and is documented in more detail in Chapter 4 of the UNCCD Good Practice Guidelines (Mattina et al., 2018). It is also implemented in the Trends Earth Plugin tool (Trends. Earth, 2019). We restrict ourselves here to the Tier 1 approach, which is based on readily available global/ regional earth observation data and geospatial modelling techniques. Starting from a baseline SOC stock map, the UNCCD-modified IPCC Tier 1 approach (IPCC, 2006) considers three change factors: 1 a land-use factor (FLU) that reflects carbon stock changes associated with land-use change; 2 a management factor (FMG) representing the main management practice specific to the land-use sector (e.g., different tillage practices in croplands); and 3 an input factor (FI) representing different levels of carbon input to soil.
Assuming land-cover can be a stand-in for land-use, the FLU change factor can be populated from land-cover and its annual transitions. The European Space Agency (ESA) CCI-LC 300 m dataset was selected as default Tier 1 data for the assessment of the land-cover trend. Because there are no known reliable global or Argentinian data at F I G U R E 2 Exponential decay function used to weigh dynamic covariates over previous years. The horizontal dashed line shows the threshold value below which exponential functions were truncated sufficient spatial resolution to obtain information on FMG and FI, these change factors were not included in the analysis.

| Model selection
Recursive feature elimination indicated that the optimal number of static covariates to be included was around 20 (Figure 3). Including more covariates did not lead to a further decrease of the RMSE. Climate covariates were most important, with precipitation covariate ranks between 1 and 12 and temperature covariate ranks between 8 and 35 (Figure 4). Elevation (rank 9) and NDVI (ranks between 13 and 33) also contributed. NDVI covariates were important for all three exponential decay parameter values, indicating that considerable temporal smoothing of NDVI was applied before it was used to predict SOC stock (see Figure 2).

| Cross-validation
A density scatter plot of predicted against observed SOC stock ( Figure 5) shows that predictions were not biased (ME was small, 0.03 kg C m −2 ) and that prediction errors were in some cases high (RMSE, 2.04 kg C m −2 ). Note that SOC stock values in the low range prevailed and that in this range, predictions and observations were concentrated along the 1:1 line. We also computed the modelling efficiency (Janssen & Heuberger, 1995). This showed that the QRF model only explained 45% of the 0-30 cm SOC stock variation. Figure 5 also shows that the predictions had less variation than the observations. This is a common characteristic of statistical prediction methods.

| Prediction maps
The spatial distribution of the 0-30 cm SOC stock was quite stable over the 1982 to 2017 time period (Figure 6). An animated GIF of the time series of annual predictions is provided in Supporting Information Figure S5. From F I G U R E 3 Root mean squared error (RMSE) for different numbers of static covariates included in the quantile regression forest (QRF) model as obtained with recursive feature elimination. Filled disc refers to the optimal number of static covariates F I G U R E 4 Variable importance for the 0-30 cm soil organic carbon (SOC) stock model (acronyms defined in Supporting Information, Table S1) [Color figure can be viewed at wileyonlinelibrary.com] Figure 6 it is difficult to detect temporal trends in the SOC stock (but see Figures 8 and 9 below). The spatial pattern largely agreed with the 2017 Baseline Neutrality Land Degradation (BNLD) 0-30 cm SOC stock map of Argentina (Supporting Information Figure S6). The largest SOC stock is observed in the mountainous regions along the border with Chile and in the very south of Patagonia. The lowest SOC stock is found in other parts of Patagonia, whereas higher SOC stocks occur in the western part of the country, such as in the Pampa, southwest of Buenos Aires.
The RMSE statistic shown in Figure 5 indicated that prediction uncertainties were large. This was confirmed by the quantile maps (shown in Figure 7 for the year 2017; results for the other years are similar). The difference between the 0.05-and 0.95-quantile SOC stock maps is large, indicating that the 90% prediction interval was wide for the whole country.

| Time series of spatial aggregates
For the entire country of Argentina, the average 0-30 cm SOC stock showed a small downward trend over time, although it stabilized between 1990 and 2000, and slightly increased after 2010 (Figure 8). The overall decrease from 1982 to 2017 is 0.0758 kg C m −2 , which equals 210.7 Gg (3.0% of the total 0-30 cm SOC stock in Argentina). Time series of the 0-30 cm SOC stock were also derived for three eco-regions of Argentina, whose geographic extent is shown in Supporting Information Figure S7. The Pampa region presented the highest mean SOC stock and had an overall decreasing trend, with a temporary increase from 1995 to 2002. The Espinal region had the lowest mean SOC stock of all three eco-regions considered. Its SOC stock gradually decreased from 2.86 kg C m −2 in 1982 to 2.69 kg C m −2 in 2017. The Chaco Seco region had stable SOC stocks over time.

| Comparison with UNCCD-modified IPCC Tier 1 results
We also compared the predicted 0-30 cm SOC stock change map between 2001 and 2015 with that obtained using the UNCCD-modified IPCC Tier 1 method (Figure 9).
• According to the UNCCD Tier 1 method, there was hardly any change in the topsoil SOC stock. • Based on the machine learning model, there are areas where the topsoil SOC stock changed substantially, with in some cases predicted changes (both up and down) of more than 1.2 kg C m −2 . • According to the machine-learning approach, the topsoil SOC stock for the entire country, over the 15-year period, decreased by 153.4 Gg C, which is seven-fold larger than the 22.7 Gg C estimated by the UNCCD Tier 1 method.

| DISCUSSION
4.1 | Space-time modelling has high uncertainty using currently available data The spatial patterns of SOC changes obtained by the QRF machine learning model agreed with recently published baseline maps of Argentina developed using "conventional" spatial DSM (Guevara et al., 2018), although these baseline maps had systematically higher SOC stocks. The systematic difference may be caused by the fact that we predicted the median SOC stock, whereas these studies derived the mean SOC stock (which is greater than the median for right-tailed distributions). Another cause of differences is that we excluded SOC data from before 1982 for model calibration, which is a substantial subset of the total dataset (Supporting Information Figure S2). All of the tested dynamic covariates to account for the delayed effect of vegetation changes on SOC contributed to the model performance. Performance may further improve if more flexibility is applied and different decay functions are used. In particular, the model could be extended with other dynamic covariates, such as dynamic climate variables.
The added value of the approach taken in this paper compared to static DSM studies is that in addition to spatial patterns, we also derived temporal trends in SOC stock, with annual prediction maps for a 36-year time period. These trends appear realistic, but should be interpreted with care. The predicted changes in the 0-30 cm SOC stock over time at point locations were an order of F I G U R E 6 Prediction maps of the 0-30 cm soil organic carbon (SOC) stock (kg C m −2 ) for the years 1982-2017, here shown in 5-year time increments (all 36 maps are shown as an animation in the Supporting Information) [Color figure can be viewed at wileyonlinelibrary.com] magnitude smaller than the uncertainty associated with these predictions. Apparently, the signal-to-noise ratio was too weak to get sufficiently narrow prediction intervals. This may not come as a surprise, in view of the scales considered here. SOC temporal variation is much smaller than SOC spatial variation, by over one order of magnitude. Being based on a calibration set of "only" 5,114 SOC stock observations, computed for points (profiles) unevenly sampled in space and time, one cannot expect highly accurate results at point locations. Both quantile regression forest and 10-fold cross-validation showed substantial uncertainty, where the latter may F I G U R E 6 (Continued) have underestimated the uncertainty somewhat. Default cross-validation tends to yield too optimistic results in the case of clustered sampling, in which case spatial or spatiotemporal cross-validation may provide more realistic uncertainty estimates (Meyer, Reudenbach, Hengl, Katurji, & Nauss, 2018;Roberts et al., 2017). Note that additional uncertainty was introduced because bulk density measurements were available for only 10% of the sites sampled for SOC. The PTF for bulk density produced biased predictions and had a large random error component, although part of the differences shown in Figure 1 are likely caused by measurement errors in the bulk density validation data. The tendency to overestimate bulk density may have led to a systematic overestimation of the Argentina 0-30 cm SOC stock.

| Averaging predicted SOC stocks over regions and the entire country
Although the uncertainty at point locations is large, predictions become less uncertain when averaged over regions. This is because within a given region both positive and negative prediction errors occur, which are cancelled out when spatial averages are taken (Heuvelink, 1998, Section 2.5). This is why we computed time series of spatial averages for the whole country and three selected eco-regions. Although not quantifiable without modelling the spatiotemporal correlation of prediction errors, the trends shown in Figure 8 have lower uncertainty than the point support predictions shown in Figure 7. Although one must take care with interpreting these time series, some tentative observations can be made. The overall decreasing trend in SOC stock for the Pampa region over the 36-year time period may be explained by the loss of natural grasslands and pastures during this period (Viglizzo et al., 2011). After dropping from 1982 to 1993, the time series showed an increase from 1993 to 2001, followed by a drop of the same magnitude until 2009, after which it is stable until 2017. These changes may have multiple causes, for example the implementation of zero-tillage in cropping systems, a widespread practice from 1993 onwards (Nocelli Pac, 2017), and the expansion of the soybean cropping area at the expense of grasslands (Piquer- Rodríguez et al., 2018).
The Espinal region showed a gradual decrease of the SOC stock in the period considered. According to Viglizzo et al. (2011), this region has been dominated by grasslands, which have been slightly reduced at the expense of annual crops. Also, Piquer- Rodríguez et al. (2018) showed that most of the land-use change in this region between 2000 and 2010 was from grazing land to cropland.
The stable SOC stock level between 1982 and 2017 in the Chaco Seco region seems to be not in concordance with the land-use change; in fact, during the last decades, this region has been listed as a global deforestation hotspot. According to the Argentinian Ministry of Environment and Sustainable Development, the annual deforestation rate in the Chaco region varied between 170,000 and 320,000 hectares between 2001 and 2013. Baldassini and Paruelo (2020) showed that land-use change in this region has generally had a negative effect on SOC stocks, but not in all cases. It also depends on land management. They also showed that conversion from shrublands to grazing areas might even increase SOC. This may explain why the predicted SOC stock in this region was stable over time, although a more thorough analysis is needed to investigate this explanation. Given the substantial uncertainty in our predictions, we cannot rule out that the Chaco Seco region did experience a SOC stock decrease during the time period considered, but that we were not able to detect it with the available data.
4.3 | Machine-learning approach is more sensitive than UNCCD Tier 1 approach The lack of observed change in SOC stock between 2001 and 2015 using the UNCCD-modified IPCC Tier 1 approach may be a result of the fact that we ignored the land management and input factors when computing these maps, because of lack of information. By ignoring land management and input factors, SOC stock change could only occur where land-use change occurred. In reality, there may be large variations in SOC stock and important changes in SOC stock over time within the same land-use type. This suggests that the SOC stock change as predicted here using the UNCCD Tier 1 approach underestimates the real change.

| Recommendations for improvement of predictions
Although spatial aggregation reduces uncertainty, from a policy and decision-making perspective it is desirable to assess SOC stock change at high spatial resolution. For instance, this is crucial to evaluate the effectiveness of carbon sequestration projects in specific regions, and to identify regions within countries and agro-ecological zones where unexpected increases and decreases of SOC stocks occur (e.g., "hot spots"), and relate this to land-use and land-management change. Sufficiently accurate information at high spatial resolution may be obtained using local SOC monitoring campaigns (e.g., De Gruijter et al., 2016), but ideally a more global coverage of SOC stock change at high spatial resolution is obtained, such as using the approach presented in this work. For this, it is imperative that the prediction accuracy is improved.
This study was an initial attempt to extend DSM from spatial to space-time SOC prediction. We identify four main ways to substantially improve prediction accuracy.
1 Collect more SOC observations. This is an obvious recommendation but that does not make it less true. One possibility is to enlarge the study area and incorporate high-quality SOC data from other countries. Note that machine learning models can benefit from calibration data that are collected in very different parts of the world, provided the environmental conditions (i.e., the soil forming factors) are similar. Thus, prediction accuracy is likely to improve if we add soil data from other parts of the world that have similar conditions to the data-poor areas (e.g., the Andes). The Argentina dataset was also "deficient" in the sense that it contains no repeated measurements over time at the same locations (e.g., data from monitoring networks), which still is the case for many countries (Smith et al., 2020). Modelling temporal change can benefit greatly from such repeated measurements, because it filters out spatial variation. In recent years, various countries and organizations have implemented soil monitoring campaigns (e.g., Gubler, Wächter, Schwab, Müller, & Keller, 2019;Orgiazzi, Ballabio, Panagos, Jones, & Fernández-Ugalde, 2018;Smith et al., 2020); making these data available for modelling could greatly improve prediction accuracy. Such monitoring schemes are becoming much easier to practically implement with the use of soil sensing technologies, as reviewed in Smith et al. (2020). However, it should be noted that for SOC stock mapping we not only need measurements of SOC concentration, but also bulk density and proportion of coarse fragments. Apart from repeated measurements over time, the spatial distribution of the sampling locations may also be optimized for random forest (Wadoux, Brus, & Heuvelink, 2019). 2 Include more relevant covariates. The accuracy of machine learning DSM models is largely determined by the predictive power of the covariates. In the case of space-time DSM, it is essential that these covariates capture the dynamics of the soil properties considered. In this study we only used a vegetation index as dynamic covariate, albeit for multiple seasons within a year and using different exponential decay functions.
In particular, climate-related covariates could potentially improve prediction, but time series on land management or other land-use-related variables are important too. In our case, we limited this study to dynamic covariates that are globally available (recall that we used the Argentina case as a pilot to test a method to be applied globally). Starting from around the year 2000 instead of 1982 would allow inclusion of many more remote sensing sources, whereas regional applications could benefit from dynamic covariates specifically available for those regions, such as localscale agricultural information on cropping area and crop production. 3 Use better models. Within the machine learning framework, we might explore recent model improvement approaches such as deep learning (e.g.,  and random forest spatial interpolation (Sekuli c, Kilibarda, Heuvelink, Nikoli c, & Bajat, 2020). Regression kriging (i.e., kriging of the residuals of the machine learning model) will also improve prediction accuracy (Hengl, Heuvelink, & Stein, 2004), at least in areas and time periods that have higher sampling density. Note that this would require the use of space-time kriging, the theory of which and software packages for are well-developed (Gräler, Pebesma, & Heuvelink, 2016). Mechanistic modelling approaches (see the Introduction and Paustian et al., 2019) potentially can get around some of the limitations of the space-time DSM approach. However, they also have disadvantages, such as lack of detailed information about model inputs, and initial and boundary conditions, but also here remote sensing might provide solutions and generate the required information. Hybrid approaches, such as machine learning on the residuals of mechanistic models, might be a valuable approach too. 4 Aggregate predictions over space and/or time. We already indicated that averaging predictions over space and/or time reduces uncertainty because positive and negative errors within a region and/or time period are cancelled out. Thus, if we are willing to sacrifice resolution in space and time, we can more quickly reach required accuracy levels. Although intuitively clear, this can also be shown mathematically using geostatistical approaches. The theory of block kriging is well developed and quantifies the reduction of uncertainty caused by spatial aggregation (Webster & Oliver, 2007 Section 4.8). In our case, this would require that the space-time semivariogram of the residuals of the machine learning SOC stock model is known, after which space-time regression block kriging could be applied.
The proposed solutions for reducing prediction uncertainty are computationally challenging and major investments are needed to narrow down prediction uncertainty to fit specific "user-demands". In particular, standardized soil sampling and soil monitoring schemes are needed (De Gruijter et al., 2016;Gubler et al., 2019;Orgiazzi et al., 2018;Smith et al., 2020), although these only offer solutions for the future, not for the past and present. If in addition the resulting data are shared in a common database (Paustian et al., 2019) and covariate data are easily available, then jointly we should be able to get a reliable system in place that informs on the status and trends in soil organic carbon and other key soil properties. Our ambition is to ultimately visualize SOC changes in near real time, as is currently possible for monitoring forest biomass change (Baccini et al., 2017). This is a true challenge, because forests are much more easily observed through remote sensing technology than soils, but given sufficient resources it is feasible. This could be embedded in a strategic research agenda for monitoring SOC change (Bray et al., 2019).

| CONCLUSION
In this paper we extended machine learning for digital soil mapping from spatial to space-time applications. From a methodological point of view this is not a complicated task because the algorithms used do not change fundamentally, but it is imperative that all soil calibration data have a time stamp as well as geographic coordinates, that dynamic covariate maps that have sufficient predictive power are available for the entire prediction period, and that computational resources are adequate to train the model and run it for all space-time prediction points. Further, a dense, well-distributed coverage of soil data in space, time and covariate space is important.
This research showed that major and sustained investments in soil monitoring, database development, covariate selection and modelling are needed to reduce prediction uncertainties and detect statistically significant SOC stock changes at high spatial resolution. Many countries and organizations across the world have implemented, or are implementing, soil monitoring campaigns. Prediction accuracy of space-time SOC stock models will dramatically increase if this becomes the norm and the data are shared for modelling.