Spatiotemporal dynamics and drivers of bare soil albedo in European croplands

Bare soil albedo plays a critical role in regulating surface energy balance and land-atmosphere interactions in agricultural systems, yet its spatiotemporal variability and controlling factors remain poorly quantified at the field scale across heterogeneous cropland landscapes. To address this, we investigate the spatial patterns and temporal dynamics of bare soil albedo in European croplands. We develop a method to reconstruct field-scale, spatiotemporally continuous bare soil albedo at a 5 d temporal resolution and 0.3 km spatial resolution using Sentinel-2 reflectance observations. Bare soil periods are identified by multiple spectral indices and the corresponding soil albedo values are derived for the period 2018–2020 using a novel machine learning framework. Two random forest models were employed to separately capture the long-term spatial structure and short-term temporal anomalies of bare soil albedo, allowing gaps caused by clouds, snow, and vegetation cover to be bridged. Model evaluation against independent site observations and existing products shows that the estimated bare soil albedo reproduces observed spatial gradients and seasonal variability across European croplands. Such variations in bare soil albedo are jointly controlled by soil properties, observation geometry and short-term soil moisture dynamics rather than by any single factor. Because these variations are of the same order of magnitude as radiation management solutions, soil radiative properties must be considered in their assessment. The resulting bare soil albedo offers a process-oriented basis for improving the representation of surface radiative properties and land-atmosphere coupling in agroecosystem and land surface models.

1 Introduction

Albedo, the proportion of incoming solar radiation reflected by a surface, modulates the local and global heat budgets by controlling the amount of incoming solar radiation that is absorbed and transformed into heat. Given that croplands occupy 37 % of the Earth's land area (FAO, 2023), it is important to quantify the albedo-related climate impacts of agricultural practices (Davin et al., 2014; Sieber et al., 2022) as they can offset or enhance the effects of greenhouse gas emissions (Bright et al., 2015).

In agricultural areas, the soil is bare (i.e., without any plant or litter cover) between successive growing periods, typically during a few months between the harvest of one crop and the sowing of the next crops. Bare soil albedo affects crop growth by modulating the local surface energy balance (e.g., through soil temperature and soil moisture) (Bsaibes et al., 2009; Zhang et al., 2022). It also influences the magnitude of agriculture-induced albedo changes and the resulting climate impacts (i.e., radiative forcing). Bare soil albedo is controlled by a combination of intrinsic soil properties and varying environmental conditions (He et al., 2019). Soil properties, including texture, mineral composition, roughness, organic matter content, and bulk density determine the optical characteristics of soil surface (Matthias et al., 2000; Cierniewski et al., 2018). As a result, they affect soil albedo directly and indirectly (via modulating the relationship between albedo and water availability). Soil moisture and solar illumination geometry further modulate the short-term temporal variability of soil albedo by altering the amount and direction of the observed reflected solar radiation (Hapke, 1981; Gascoin et al., 2009; Badía et al., 2016). The interaction of these factors, as well as the frequent agricultural activities (e.g., harvest, changes of crop types), gives rise to pronounced spatial heterogeneity and rapid temporal variability, particularly in intensively managed agricultural landscapes (Starr et al., 2020; Sieber et al., 2022). Nevertheless, the spatiotemporal variability and controlling mechanisms of bare soil albedo remain insufficiently quantified at field scale across heterogeneous cropland regions.

Ground-based measurements, airborne or unmanned aerial vehicle measurements measure incoming and outgoing solar radiation over the soil surface, but have limited spatial coverage (Canisius et al., 2019; Pastorello et al., 2020; Amazirh et al., 2021; Pan et al., 2024). Satellites observations provide comprehensive spectral reflectance data with extensive temporal and spatial coverage, allowing for the monitoring of bare soil albedo datasets at regional to global scale (Lesaignoux et al., 2009; Verheijen et al., 2013; Carrer et al., 2014; Cierniewski et al., 2017; Rizzo et al., 2023). However, commonly used satellite-based products, such as those derived from MODIS, typically operate at spatial resolutions of several hundred meters to kilometers and temporal resolutions of weeks to months, which limits their ability to resolve field-scale heterogeneity and short-term albedo dynamics driven by soil moisture fluctuations and agricultural management (Zhou et al., 2005; Houldcroft et al., 2009; Cheng and Liang, 2014; He et al., 2019).

For example, the ECOCLIMAP project provides bare soil albedo datasets at a coarse spatial resolution of 1 km and a temporal resolution ranging from 10 d to monthly, which is designed for use in climate modelling, but does not meet the resolution requirements for effective agricultural monitoring able to capture the impacts of agricultural management (Masson et al., 2003; Faroux et al., 2013; Carrer et al., 2014). While products at sufficiently high spatiotemporal resolution exist but are limited to small regions or short time periods, constraining their applicability for continental-scale process analysis (Bsaibes et al., 2009; Roujean et al., 2018).

Numerical models based on physical processes such as radiative transfer models are also applied to produce bare soil albedo combined with satellite-based reflectance (Jacquemoud et al., 1992; Bablet et al., 2018). For example, a global ∼50 km monthly bare soil albedo map was generated by radiative transfer equations from 2001 to 2009 based on MODIS surface albedo (Pinty et al., 2011). However, this coarse spatiotemporal resolution limits its utility in evaluating albedo-related climate impacts of land management. Additionally, because of issues related to calibration and validation, the reliability of the simulations provided by process-based models is limited (Bsaibes et al., 2009).

Machine learning can investigate bare soil albedo dynamics by considering complex dependencies between variables and deliver more accurate and reliable predictions than traditional statistical methods (Jones, 2017; Nguyen, 2021; Buster et al., 2024). This capability is particularly valuable for understanding the complex interactions between bare soil albedo, climate, and pedological features, as opposed to traditional statistical approaches that typically rely on simple relationships between bare soil albedo and single features such as the exponential relationship between bare soil albedo and soil water content (Lesaignoux et al., 2009; Gascoin et al., 2009; Verheijen et al., 2013). Additionally, machine learning-based predictions of bare soil albedo can address spatiotemporal gaps in satellite observations (Ren et al., 2022) caused by snow and cloud cover, enabling large-scale continuous and dynamic analysis and simulation of bare soil albedo in crop regions under historical, present and future scenarios. When combined with high-resolution satellite imagery (e.g., Sentinel-2 satellite with 0.01 km spatial resolution), machine learning enables the reconstruction of continuous bare soil albedo at both large scales and high resolution, especially in cropland regions.

In this study, we investigate the spatiotemporal variability and drivers of bare soil albedo across European croplands by integrating Sentinel-2 observations with a machine-learning reconstruction framework. We test two hypotheses: (1) for the spatial distribution of bare soil albedo, soil properties exert a stronger control than radiative factors, and (2) for temporal variability of bare soil albedo, fluctuation in soil moisture is more influential than static soil properties. To address these hypotheses, our objectives are to (1) quantify the spatial patterns of bare soil albedo at the field scale across Europe, (2) characterize the temporal dynamics of bare soil albedo during bare soil periods, and (3) identify the relative contributions of soil and topographic properties versus dynamic drivers such as soil moisture and solar geometry. By doing so, we provide a process-oriented assessment of bare soil albedo variability that is relevant for understanding the energy balance of agricultural systems and improving the representation of surface radiative properties in land surface models.

2 Material and methods

2.1 Data source and processing

Our approach consists of the following steps (Fig. 1): First, we identified bare soil periods and extracted corresponding bare soil albedo (α_bs) values from Sentinel-2 surface albedo products based on thresholds of satellite-derived indices. Second, we trained two random forest (RF) models on Sentinel-2 surface albedo during bare soil periods for the prediction of spatial distribution and temporal anomaly of α_bs at the European cropland scale. Third, we compared the predicted α_bs to existing α_bs datasets at site or European scales.

Figure 1 The flow chart of the data processing procedure. “NDVI, NBR2, GI, NBR” are the normalized difference vegetation index (NDVI), normalized burn ratio 2 index (NBR2), greenness index (GI) and the normalized burn ratio index (NBR), respectively. “LSZA, SM, SM_RC, SOC, bulkdens” represent local solar zenith angle (LSZA), soil moisture (SM), relative SM anomaly (SM_RC), soil organic carbon (SOC) and bulk density (bulkdens), respectively. “α_bs,spatial”, “α_bs,temporal” and “α_bs,RF” are the spatial distribution and temporal variation of bare soil albedo (α_bs) predicted by random forest models, and the final spatiotemporal α_bs maps, respectively. R² and RMSE are the coefficient of determination and root mean square error. The bold black dotted (“1”) and solid arrow (“2”) lines represent the procedure of spatial and temporal α_bs prediction, respectively.

2.1.1 Daily surface albedo

We derived 5 d 0.3 km black-sky albedo (BSA) and white-sky albedo (WSA) from the 0.01 km Sentinel-2 surface reflectance product for the period 2018–2020 using the approach of Lin et al. (2022). This approach uses a look-up table trained on prior information on the angular-based relationships between surface reflectance and broadband albedo derived from the 0.5 km resolution MODIS BRDF product (MCD43A1). The look-up table is used to derive BSA and WSA from Sentinel-2 surface reflectance product assuming the same relationship as between BSA and WSA and MODIS BRDF product. The daily albedo derived from Sentinel-2 reflectance shows a higher correlation with site observations (n=588) (R²=0.44, RMSE = 0.004) compared to the MODIS product (R²=0.16, RMSE = 0.06) (Fig. S7 in the Supplement). This provides support for the use of the Sentinel-2 reflectance in combination with the look-up table trained on MODIS.

A cloud probability product was utilized to filter out images affected by cloud cover. This process employed cloud masking parameters as outlined in the Google Earth Engine Sentinel-2 cloud masking tutorial (https://developers.google.com/earth-engine/tutorials/community/sentinel-2-s2cloudless, last access: 31 October 2026).

Specifically, the “CLOUD FILTER” and the “BUFFER” parameters were set to 10 % and 0.02 km in this process. This framework was accomplished on the Google Earth Engine platform to produce the datasets of BSA and WSA with a spatial resolution of 0.3 km at a European scale. We used 0.3 km as the spatial resolution to match the average field size in Europe and ensure a sufficiently large number of data points required for machine learning (see below).

We finally obtained the 5 d shortwave surface albedo (α) from BSA and WSA weighted by the fraction of diffuse solar radiation (k_d):

$$\begin{array}{}\text{(1)}& \mathit{\alpha }=\left(\mathrm{1}-k_{\mathrm{d}}\right)\times \text{BSA}+k_{\mathrm{d}}\times \text{WSA}\end{array}$$

The parameter k_d is computed based on the global surface solar radiation (S₀, MJ m⁻² d⁻¹) and extra-terrestrial solar radiation (S_e, MJ m⁻² d⁻¹) during the same periods with BSA and WSA datasets. Details of the diffuse solar radiation model can be found in Yang et al. (2019). S₀ corresponds to the daily average of hourly downward surface solar radiation at the spatial resolution of ∼11.1 km from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis 5 (ERA5) (Hersbach et al., 2020). S_e is estimated from the solar constant, the solar declination and the time of the year in a specific location, which is described by Allen et al. (2006).

Filtering Sentinel-2 surface albedo during bare soil periods from 2018 to 2020 using spectral indices, we obtained 132 487 soil samples in the European cropland regions. The average Sentinel-2 α_bs of all 0.3 km resolution bare soil pixels is approximately 0.18±0.03 (mean ± standard deviation). Details of extracting bare soil periods in European cropland regions and defining the 0.3 km spatial resolution of Sentinel-2 α_bs are described below.

We computed (1) the temporal average at each 0.3 km sample location for training and testing an RF model for the spatial pattern of α_bs, and (2) the temporal anomaly of α_bs as the difference between the original α_bs and the temporal average. Samples with an anomaly of 0 were removed as they represent locations with only a single α_bs observation, thus lacking temporal representativeness. The α_bs anomaly dataset was used as the dependent feature for training and testing RF for the temporal anomaly of α_bs.

2.1.2 European cropland mask

The ESA WorldCover product at 0.01 km spatial resolution generated from Sentinel-1 and Sentinel-2 observations in 2021 was utilized to identify croplands in Europe (https://esa-worldcover.org/, last access: 31 October 2026) (Zanaga et al., 2022). This time-invariant product includes 11 land cover types, including one class for croplands.

2.1.3 Average crop field size

We leveraged the Land Use/Cover Area frame Survey (LUCAS) 2018 Copernicus land cover dataset (d'Andrimont et al., 2021) to estimate the averaged field size in European croplands to decide on the spatial resolution. The LUCAS dataset encompasses 63 287 polygons distributed across various land cover types with 120 attributes in Europe, including artificial lands (0.6 %), croplands (34.1 %), woodlands (30.0 %), shrublands (4.5 %), grasslands (28.0 %), bare lands (2.5 %), wetlands (0.3 %) and water (0.0 %). The field sizes of cropland in the LUCAS dataset were labelled as categorical attributes in four levels, including “>10 ha”, “1–10 ha”, “0.5–1 ha” and “<0.5 ha”. We assigned representative values of 15, 5, 0.75 and 0.25 ha for each category (Fig. S8) to approximate the averaged radius of external circles of cropland regions at 0.16 km, varying from 0.03 to 0.22 km (n=36 199). For simplification, 0.15 km radius (pixel size of 0.3 km) was finally defined as the spatial scale of the Sentinel-2 α_bs dataset in this study. Aggregation of 0.01 km observations to 0.3 km is also a strategy to increase temporal Sentinel-2 coverage. While it loses fine-scale details (like hedges, within field heterogeneity), it not only increases the likelihood for suited conditions for bare soil observations at a given location but also gains robustness, reduces noise, and facilitates processing.

2.1.4 Periods of bare soil exposure

To extract periods with bare soil in European cropland regions, we used four spectral indices and their thresholds from previous analysis (Rizzo et al., 2023), including greenness index (GI) (GI <0.65), normalized difference vegetation index (NDVI) ($-\mathrm{0.05}<\text{NDVI}<\mathrm{0.25}$), normalized burn ratio index (NBR) (NBR $>-\mathrm{0.23}$) and normalized burn ratio 2 index (NBR2) ($-\mathrm{0.05}<\text{NBR2}<\mathrm{0.15}$). The thresholds of these four indices were determined based on the spectra of 8005 soil samples collected globally (Rizzo et al., 2023). In this analysis, these indices were computed from the Sentinel-2 reflectance product on GEE (Table 1). We further utilized the NDVI and normalized difference snow index (NDSI) to eliminate the periods influenced by water or snow cover on soils. Dates with NDVI lower than 0 were excluded to avoid water bodies, rocks or other non-vegetated land surfaces (Korchagina et al., 2020). An NDSI threshold of 0.1 was utilized to obtain snow-free bare soil periods, which is stricter than the commonly used threshold of 0.4 (Riggs et al., 2016; Härer et al., 2018). Bare soil period is identified if the conditions of all indices are satisfied. Finally, data points outside mean ± two standard deviations were considered outliers and removed.

Table 1 The vegetation index for extracting the bare soil periods in European cropland regions from 2018 to 2020.

The Blue, Green (VIR), Red, NIR, SWIR1 and SWIR2 are represented by bands 2, 3, 4, 8, 11 and 12 of Sentinel-2 products, respectively.

2.2 Random Forest models for mapping bare soil albedo

2.2.1 Rationale and design choices

The Sentinel-2 time series over Europe contains frequent gaps due to cloud and snow cover, which hinder regular sampling of bare soil conditions. To obtain spatiotemporally continuous fields, we used a RF model to learn the relationship between observed α_bs dynamics and potential drivers, thereby extending Sentinel-2 observations in time and space. RF is well-suited to nonlinear responses, mixed predictors, and heterogeneous sampling, and is widely applied in remote sensing regression and gap filling (Breiman, 2001).

2.2.2 Model inputs and outputs

Two RF models were developed. The first model aimed at predicting the spatial distribution of the temporal average of α_bs (α_bs,spatial). Predictions were performed from the following input features: temporal-averaged local solar zenith angle (LSZA), soil organic carbon (SOC), sand and clay contents, bulk density, pH, temporally averaged soil moisture (SM) and elevation (Table 2).

Table 2 Features used in the analysis, including satellite properties, atmosphere, soil and surface conditions, as well as country identity.

The second model aimed at predicting the temporal variation of α_bs, defined as ${\mathit{\alpha }}_{\mathrm{bs},\mathrm{temporal}}={\mathit{\alpha }}_{\mathrm{bs}}-{\mathit{\alpha }}_{\mathrm{bs},\mathrm{spatial}}$. The input features of this second model were LSZA, time of day, days of year, relative SM anomaly (SM_RC) and the predicted values of α_bs,spatialprovided by the first RF model (Table 2). All the input features were linearly resampled from the original resolution to 0.3 km pixel size and 5 d temporal resolution to match the pixel size of α_bs.

The selection of features was based on theoretical and empirical evidence for their potential control on α_bs. SOC, sand and clay contents, bulk density and pH affect α_bs by modulating soil color, microtopography, surface roughness and water-holding capacity (Cierniewski et al., 2014; Sadeghi et al., 2018; Elfarisna et al., 2024). For example, light-colored sand soil generally has a higher α_bs due to larger particles with smoother surfaces, as well as weaker water-holding capacity for reflecting solar radiation, compared to clay and silt (Sun et al., 2015). The elevation, extracted from the digital elevation model (DEM), affects the distribution of incoming solar radiation on the soil and, consequently α_bs (Wen et al., 2013). The LSZA, days of year and time of day affect α_bs by controlling the direction and the path length of solar radiation in the sun-surface-satellite geometry (Yang et al., 2008). SM is negatively correlated with α_bs by changing soil colors and influencing the soil roughness for specific soil types (e.g., dry soil tends to form clods or aggregates, increasing surface roughness) (Gascoin et al., 2009; Sugathan et al., 2014; Rizzo et al., 2023).

The five 0.25 km time-invariant top (upper 30 cm) soil properties (SOC, sand and clay contents, bulk density and pH) were extracted from the U.S. Department of Agriculture (USDA) soil datasets of OpenLandMap (Hengl, 2018a, b, c; Hengl and Wheeler, 2018). The 0.03 km static DEM was obtained from the Shuttle Radar Topography Mission (SRTM) (Farr et al., 2007). The satellite properties (LSZA, days of year and time of day) were derived from Sentinel-2 reflectance products with a 5 d interval and 0.3 km spatial resolution. The SM dataset is from the 1 km daily Global Surface Soil Moisture (GSSM) product (Han et al., 2023), which provides the temporal dynamics of the daily surface soil moisture (0–5 cm) in cm³ cm⁻³, at a coarse spatial resolution of 1 km in European scale from 2000 to 2020. GSSM was predicted from machine learning using the datasets of the International Soil Moisture Network, satellite and meteorology (Han et al., 2022).

The SM_RC is calculated as the relative deviation from the temporal-averaged SM in each pixel:

$$\begin{array}{}\text{(2)}& {\text{SM\_RC}}_{i,j}={\displaystyle \frac{{\text{SM}}_{i,j}-{\text{SM}}_{\mathrm{ave},j}}{{\text{SM}}_{\mathrm{ave},j}}}\end{array}$$

Where SM_i,j and SM_RC_i,j are the SM value and relative SM anomaly at time i in pixel j, respectively. SM_ave,j is the temporal-averaged SM in pixel j.

2.2.3 Training and testing of the models

The procedure involves hyperparameter tuning, model training and testing. This procedure was implemented to build the two RF models predicting α_bs,spatial and α_bs,temporal.

We used a dataset including a total of 75 941 values of the temporal average of α_bs for training and testing the first model, and a dataset including a total of 56 754 values of the temporal anomaly of α_bs for training and testing the second model. Each dataset was randomly split into two samples: one for training and hyperparameter tuning (80 %), and one for out-of-box testing (20 %). Hyperparameter tuning was conducted for each model separately, using 10-fold cross-validation with grid searching on the training dataset. Hyperparameters include the number of decision trees (ntrees) from 300 to 1000, and the number of variables sampled (mtry) from 2 to 7 at each split. Specifically, for each pair of hyperparameters (e.g. ntrees of 300 and mtry of 2). The procedure of 10-fold cross-validation is described in Sect. S1 in the Supplement.

With the selected optimal hyperparameters, the RF model was trained on the whole training dataset and hyperparameter tuning sample, and assessed on the out-of-box testing sample using R² and RMSE as quality indicators. Overfitting was additionally checked by comparing the final R² against the R² obtained during cross-validation.

The SHAP (SHapley Additive exPlanations) values indicate the contribution of each value for each feature to the individual model prediction by considering all possible combinations of feature subsets (Zhang et al., 2022). The SHAP values of the features in spatial and temporal RF models were induced to visualize the dependencies among features in jointly deciding the α_bs variation.

2.2.4 Predicting bare soil albedo over European croplands

The trained spatial and temporal RF models were then applied to produce the 5 d 0.3 km European α_bs datasets based on the selected features: (1) LSZA, SM, elevation, SOC, pH, clay content, sand content and bulk density for predicting α_bs,spatial with the spatial RF model; (2) LSZA, α_bs,spatial, SM_RC, days of year, and time of day for predicting α_bs,temporal with the temporal RF model. The final continuous spatial-temporal predictions (α_bs,RF) were derived by adding the predictions of α_bs,temporal and α_bs,spatial, as follows

$$\begin{array}{}\text{(3)}& {\mathit{\alpha }}_{\mathrm{bs},\mathrm{RF}}={\mathit{\alpha }}_{\mathrm{bs},\mathrm{spatial}}+{\mathit{\alpha }}_{\mathrm{bs},\mathrm{temporal}}\end{array}$$

2.3 Independent evaluation datasets

Five types of α_bs datasets were utilized to evaluate the performance of the final spatiotemporal albedo predictions (α_bs,RF) at seven ICOS European cropland sites (BE-Lon, FR-Lam, FR-Aur, FR-EM2, DE-RuS, DE-Geb, DE-Kli), one European Fluxes Database Cluster (EFDC) European cropland site (CZ-KrP), and one Swedish cropland site in Uppsala, and the overall European cropland regions during bare soil periods (Fig. 3, Table 3). The bare soil periods at ICOS and EFDC sites were manually identified by checking daily or weekly site photos from 2018 to 2020, selecting periods when the soil within the visual range of the camera was visibly exposed without any covers (Fig. S2). The bare soil periods at Swedish cropland sites were obtained from records directly provided in the literature.

Table 3 Evaluation datasets of bare soil albedo (α_bs) at sites or European cropland regions from the literature.

The α_bs datasets of all sites, derived either from site eddy covariance measurements or 5 d 0.3 km Sentinel-2 surface albedo observations, are denoted as α_bs,site and α_bs,S2 to represent the α_bs we used for evaluation at site scale in the following text. The α_bs,site at the Swedish site, composed of six cropland fields, was extracted from literature directly for comparison (Sieber et al., 2022). The selected eight sites encompass a broad pedoclimatic and management gradient across European croplands (Table S1 in the Supplement). They represent contrasting soil conditions, including Luvisol, Fluvisol, Cambisol, and Chernozem/Phaeozem, with the ranges of silt, clay and sand contents from 26.6 % to 78.1 %, from 14.7 % to 53.0 %, and from 7.2 % to 38.0 %, respectively. These sites also cover diverse climatic conditions in air temperature, incoming shortwave solar radiation, precipitation and cloud cover at annual scale. Seasonal variability in solar radiation and soil moisture (characterized by the seasonal minimum and maximum values in Table S1) captures the range of environmental conditions experienced across seasons. In addition, these sites implement different agricultural management practices, including crop rotations of 2–8 years and bare soil periods occurring mainly in spring, autumn and winter (approximately 30–180 d, Fig. S2).

We aggregated the field data to represent the α_bs,site for this site. The comparison at site scale is only implemented if the sample size at each site is more than three.

We additionally examined the relationships between the daily soil sensor SM measurements and different α_bs datasets to evaluate their consistency in expressing the dependence of α_bs on SM, including α_bs,RF, α_bs,MODIS, α_bs,Gascoin, α_bs,S2 and α_bs,site (Table 3).

3 Results

3.1 Controls on spatial and temporal variability of bare soil albedo

The random forest models successfully reproduced both the spatial patterns and temporal variability of α_bs across European croplands (Fig. 2a and c), providing a basis for diagnosing the dominant controls on α_bs at different scales. The spatial model achieved an R² of 0.46, while the temporal anomaly model explained 35 % of the observed variability, with both models yielding RMSEs of 0.02 based on out-of-bag evaluation datasets.

Figure 2 The performance of the two random forest models (RF) for predicting European bare soil albedo (α_bs,RF) in time and space. (a) and (c) are the scatter plots between the predicted and observed spatial bare soil albedo, and between predicted and observed temporal albedo variations from the test dataset, respectively. (b) and (d) are the importance rankings of selected features for spatial and temporal RF modelling. α_bs,spatial is the predicted spatial distribution of bare soil albedo. “LSZA, SM, SM_RC, SOC” represent local solar zenith angle (LSZA), soil moisture (SM), relative SM anomaly (SM_RC) and soil organic carbon (SOC), respectively.

The long-term spatial distribution of α_bs (α_bs,spatial) is primarily controlled by temporally averaged solar zenith angle (LSZA), mean soil moisture conditions and elevation (Fig. 2b). In contrast, temporal variability during bare soil periods is most strongly influenced by instantaneous solar zenith angle, the underlying spatial albedo pattern, and soil moisture anomalies (Fig. 2d). The importance of these variables on α_bs has been widely explored in previous studies (Idso et al., 1975; Post et al., 2000; Yang et al., 2008; Gascoin et al., 2009; Roosjen et al., 2015; Cierniewski and Ceglarek, 2023). The remaining predictors contribute to weaker but non-negligible impacts. It is reasonable that soil properties (e.g., soil organic carbon, pH) control the spatial distribution of α_bs by adjusting soil colors via the contents of organic matter and nutrients (Sirisathitkul and Sirisathitkul, 2025) while time records adjust its temporal variation (Yu et al., 2024).

Partial dependence analysis further illustrates the physical consistency of these controls. Both LSZA and SM have negative effects on α_bs,spatial and α_bs,temporal predictions (Fig. S3a, b, d and f). The elevation increases α_bs,spatial by representing specific soil types with elevation (Fig. S3c) (Badía et al., 2016; Praeg, 2020). An impact of α_bs,spatial on α_bs,temporal is expected as it represents the location-specific climatic and environmental conditions which are otherwise not represented by the selected predictors (Fig. S3e).

Figure 3 (a) The averaged bare soil albedo (α_bs) at 300 m resolution derived from the random forest models over the period from 2018 to 2020 in the European cropland regions (α_bs,RF). Dots with colors represent the temporal average of α_bs,RF at site scale; (b) The comparison between α_bs,RF and α_bs,S2 (triangles, with α_bs,S2 extracted from the Sentinel-2 surface albedo products), and between α_bs,RF and α_bs,site (circles, with α_bs,site computed from site radiation measurements) at six site fields during the same bare soil periods. The triangles and circles in different colors correspond to the averages of α_bs,RF, α_bs,S2 and α_bs,site at different sites, and the bars represent two standard deviations. The numbers in legend brackets represent the sample size (Table 3 for details). Statistics shown in black indicate the correlation passes the significance level of 0.05 (p<0.05), while the one in gray indicates non-significance (p>0.05).

3.2 European-scale spatial pattern of bare soil albedo

The spatiotemporal average of α_bs,RF is 0.15±0.01 over European cropland from 2018 to 2020 which is within the range of α_bs found in the literature, which spans from 0.15 to 0.22 derived from 7537 European topsoil samples in 2022 (Szabó et al., 2024). Our α_bs dataset has a higher spatiotemporal average than earlier products including α_bs,Carrer of 0.12±0.03 (2001–2010) and α_orc,bgd of 0.13±0.03 (2000–2009), respectively, which have coarser spatial resolutions of 1 km and ∼ 55.5 km (Figs. 4 and S4) (Carrer et al., 2014; He et al., 2019). The distribution of the final predictions α_bs,RF was similar to the distribution of the α_bs data extracted from Sentinel-2 observations by vegetation indices within European cropland regions from 2018 to 2020 (Fig. S4). This indicates that our predictions cover the spatio-temporal variations in European croplands.

Figure 4 The spatial comparison between the new albedo product α_bs,RF and two existing bare soil albedo products in European cropland scale. (a) the averaged ∼1 km 8 d European bare soil albedo dataset from 2000 to 2010 (α_bs,Carrer); (b) the difference between the averaged 2018–2020 α_bs,RF and the averaged 2001–2010 α_bs,Carrer; (c) the averaged ∼50 km monthly bare soil albedo dataset from 2001 to 2009 (α_orc,bgd); (d) the difference between the averaged 2018–2020 α_bs,RF and the averaged 2001–2009 α_orc,bgd.

The multi-year spatial pattern exhibits broad west-east and north-south gradients (Fig. 3a) that align with major soil zones, with darker surfaces where organic matter is higher and brighter surfaces where carbonates are more prevalent, a contrast repeatedly reported for European croplands (Carrer et al., 2018; Rizzo et al., 2023). The Iberian Peninsula, southwestern and southern France, and the Mediterranean coastal regions have bright soils with α_bs of 0.18–0.20. While the soil in north-central and eastern Europe (e.g., Germany and Poland) is dark with α_bs of 0.12–0.15.

Despite general agreement in continental-scale gradients, spatial correlations between our field-scale α_bs,RF estimates and existing coarse-resolution products remain low (Pearson r=0.07–0.15; Fig. 4).

3.3 Site-scale characteristics of bare soil albedo

During bare soil windows, α_bs,RF has an average of 0.15±0.01 across six eddy-covariance cropland sites, with time series showing coherent field-scale dynamics (Figs. 3b and 6). Time series of α_bs,RF exhibit site-specific background levels with short-term fluctuations.

The daily α_bs,site and α_bs,S2 at the sites, for which bare soil albedo estimates were available, were utilized to compare with daily variations in α_bs,RF (Figs. 5 and 6). Results show that the time-averaged predicted albedo was close to observation (R²=0.84, RMSE = 0.12 for comparison with α_bs,S2), while the agreement with in-site observations is weaker (R²=0.12, RMSE = 0.11 for comparison with α_bs,site) (Fig. 3b). α_bs,RF is $+\mathrm{21.87}\pm \mathrm{27.71}$ % (mean ± standard deviation) higher than α_bs,site (0.14±0.02) for the same periods (Figs. 5 and 6). While the difference between α_bs,RF and underlying training data (α_bs,S2, 0.12±0.02) across sites is smaller ($+\mathrm{7.11}\pm \mathrm{3.55}$ %). It implies that the difference between α_bs,RF and α_bs,site is a mix of biases in the RF to capture Sentinel-2 observations and differences between α_bs,S2 from α_bs,site. We also found that the difference of the model predictions between α_bs,RF and α_bs,site varies across different soil types. BE-Lon, FR-EM2 and FR-Lam sites, belonging to alfisols (light brown to reddish-brown), have lower bias compared with site measurements (Fig. 3). While DE-Geb and DE-Kli (dark brown to black molisols) and CZ-KrP (brown to yellowish-brown inceptisols) induce higher differences from measurements.

Figure 5 The boxplots of α_bs,RF, α_bs,MODIS, α_bs,Gascoin, α_bs,S2 and α_bs,site at six cropland sites during the same bare soil periods from 2018 to 2020. The values of α_bs,RF are the bare soil albedo derived from the spatiotemporal random forest models. The values of α_bs,S2 are extracted from the Sentinel-2 surface albedo products, while the values α_bs,site are computed by site eddy covariance measurements (see Table 3). The values of α_bs,MODIS are extracted from 500 m daily MODIS surface albedo product (MCD43A3). α_bs,Gascoin are computed by the exponential function of SM using the 1 km GSSM SM product (see Table 3). The red and blue lines in each figure represent the median and mean values, respectively.

Figure 6 Values of α_bs,RF, α_bs,MODIS, α_bs,Gascoin at four sites during the same bare soil periods from 2018 to 2020, compared with α_bs,S2 and α_bs,site (see Table 3). Values of α_bs,RF are the bare soil albedo derived from the spatiotemporal random forest models (red hollow dots). Values of α_bs,S2 are extracted from the Sentinel-2 surface albedo products (black hollow dots), while α_bs,site values are computed by site radiation measurements (gray hollow dots). The values of α_bs,MODIS are extracted from 500 m daily MODIS surface albedo product (MCD43A3) (yellow hollow dots). The values of α_bs,Gascoin are computed by the function of soil moisture (SM) using the 1 km GSSM SM product (brown hollow dots). “Period X” in each plot indicates the bare soil period identified by daily or weekly site photos. Plot in-between (a) and (b) belongs to BE-Lon site but in Period 2.

For reference, α_bs,Gascoin (0.17±0.01) and α_bs,MODIS (0.16±0.02) are higher than α_bs,RF and α_bs,site at the same sites and times (Fig. 5). Correlation analyses showed that temporal variations in α_bs,RF, as prescribed by Eq. (3), tracked soil moisture dynamics (R²=0.51, RMSE = 0.006), which is darkening after wetting and progressive brightening during dry-downs (Fig. 7). This correspondence is stronger than α_bs,site (R²=0.30, RMSE = 0.02) and α_bs,S2 (R²=0.06, RMSE = 0.01), comparable to α_bs,MODIS (R²=0.46, RMSE = 0.007), and weaker than α_bs,Gascoin (R²=0.94, RMSE = 0.001). The particularly high correlation for α_bs,Gascoin is expected because α_bs,Gascoin is directly derived from soil moisture. Compared with α_bs,site, α_bs,RF is consistently brighter across the soil moisture range (∼ 0.01–0.02 higher in α_bs), with the positive bias most pronounced for soil moisture >50 %. α_bs,RF also exhibits a smaller amplitude and variance than α_bs,site, smoothing the peak darkening seen in the site observations.

Figure 7 The linear relationship between α_bs,site and soil moisture (SM), α_bs,RF and SM, α_bs,MODIS and SM, and α_bs,Gascoin and SM, respectively (see Table 3). The superficial SM dataset in the x axis is derived from daily soil sensor measurements which have an SM step of 5 % during the same bare soil periods from 2018 to 2020. The dots with areas in different colors show the averages and 2-time standard deviations of all datasets in each SM step. The α_bs,RF is produced by the trained random forest models using daily GSSM SM product as one of the predictors (red dots). The α_bs,S2 is extracted from the Sentinel-2 surface albedo products (black dots), while α_bs,site is computed by site eddy covariance measurements (gray dots). The α_bs,MODIS is extracted from 500 m daily MODIS surface albedo product (MCD43A3) (yellow hollow dots). α_bs,Gascoin is computed by the function of SM using the 1 km GSSM SM product (brown dots). “*” means the linear relationship (solid lines) passes the significance level of 0.05.

4 Discussion

4.1 Characteristics of the reconstructed bare soil albedo

The reconstructed α_bs is designed to capture the dominant physical controls on bare soil albedo while balancing the trade-offs behind continuity and representativeness based on Sentinel-2 reflectance. It enables the characterization of bare soil albedo variability at spatial and temporal scales relevant for agricultural surface energy exchange across Europe. Compared with earlier continental-scale estimates (α_orc,bgd and α_bs,Carrer), the reconstructed α_bs resolves sub-kilometer heterogeneity that is critical for diagnosing field-scale radiative contrasts in fragmented cropland landscapes (Pinty et al., 2011; Carrer et al., 2014; He et al., 2019).

The differences among estimates could be attributed to differences in the periods covered: 2001–2009 for α_orc,bgd, 2001–2010 for α_bs,Carrer, and 2018–2020 for α_bs,RF. Changes in climate and environmental conditions during these periods, such as the observed decreasing trend of soil moisture on a European scale from 1991 to 2020 compared to 1991, especially in central and eastern Europe (Almendra-Martín et al., 2022), can theoretically result in average higher 2018–2020 α_bs compared to the period before 2010.

The estimation of α_bs largely reflects the seasonal variability of α_bs in European cropland fields (Fig. 6), providing a basis to quantify measures for radiation management in agriculture (Yu et al., 2024). This dataset reduces sampling discontinuities associated with cloud cover, snow, and vegetation masking (Figs. 6, S1 and S2), allowing seasonal envelopes and inter-site contrasts in α_bs to be examined more consistently (Wang et al., 2025). The selected features for RF models enable α_bs prediction in regions with limited satellite observations due to cloud cover or aerosols, because the necessary predictors of the RF model, including soil properties and solar-surface-satellite geometry, are not affected by atmospheric conditions (Fig. 2b and d).

4.2 Physical controls of the reconstructed bare soil albedo

The decomposition of bare soil albedo into spatial and temporal components provides a framework for interpreting its physical controls. Intrinsic soil properties affect bare soil albedo by establishing the optical characteristics of the soil surface and regulating soil-water interaction, while soil moisture dynamics and radiative geometry control the temporal fluctuations of the observed bare soil albedo.

Our analysis supports this framework by isolating and quantifying the contributions of soils, weather and management on spatio-temporal variations in bare soil albedo. It reveals that the α_bs dynamics arise from the interplay of multiple controls rather than any single dominant driver (Yang et al., 2008). Temporal changes (α_bs,temporal) are primarily associated with variations in soil moisture and local solar zenith angle, with additional seasonality represented by time-of-day and day-of-year terms (Fig. 3b). This underscores the role of short-term variations in illumination geometry and surface wetness in shaping bare soil albedo dynamics (Yu et al., 2024). Whereas spatial contrasts (α_bs,spatial) reflect long-term patterns associated with soil moisture and solar geometry, further modulated by invariant factors such as elevation and soil type (Fig. 3d). This suggests that stable radiative conditions and soil-related properties (e.g., soil organic carbon and pH) exert dominant controls on large-scale albedo variability (Sirisathitkul and Sirisathitkul, 2025). Approaches that rely solely on soil moisture as drivers of α_bs (Gascoin et al., 2009) are thus likely to misestimate α_bs such as overestimating α_bs during dry-down over dark, carbon-rich Mollisols, and underestimating α_bs during wet periods under low solar zenith angles, when viewing geometry elevates apparent reflectance.

Soil moisture acts as a short-term control on α_bs by changes in soil color and soil roughness through wetting, darkening and dry-down brightening (Rizzo et al., 2023). The magnitude and duration of these responses depend on the soil water retention capacity (Post et al., 2000), resulting in different α_bs responses to the same soil moisture change for different soil textures. The impact of local solar zenith angle is expected due to its impact on the relative position between the spaceborne sensor and the observed areas, affecting the amount of reflected solar radiation (Yang et al., 2008): a small zenith angle (e.g., 25°) results in nearly vertical solar radiation reflectance and a high α_bs (Fig. S3a). It is treated as a physical control on the measured radiative signal (Liu et al., 2009). Its contribution is separate from soil moisture and soil properties and reflects the presence of sun-surface-satellite geometry rather than a proxy of unavailable factors, meaning that α_bs can be understood in light of seasonal changes in observation geometry as well as surface state.

Soil properties explain much of the stable spatial structure onto which moisture and geometry act (Rizzo et al., 2023). For example, soil texture influences α_bs through both direct and indirect pathways. Particle size, mineral composition directly influences soil optical properties (e.g., surface color and scattering characteristics), thereby establishing the inherent α_bs. In addition, these properties related with texture regulate the soil hydraulic behaviour (Sadeghi et al., 2018). Variables such as elevation or bulk density do not change α_bs directly but act as proxies for topographic setting or soil structure where more specific drivers (e.g., aspect or slope at the same scale) are not available in the RF model (Fig. 2b and d). For instance, higher elevations generally have lighter soil colors, which might increase α_bs,spatial (Fig. S3c) (Badía et al., 2016; Praeg, 2020). Because of the collinearity of bulk density and clay content (Heuscher et al., 2005), their importance reflects predictive utility rather than causal influence on α_bs. They stabilize spatial patterns and improve field coherence, not to assign mechanistic primacy. A day-of-year term is included to capture recurring within-year patterns that arise from sunlight, crop calendars, and hydro-meteorology (Reichstein et al., 2019). Since the same day-of-year repeats across years, this term is not intended to encode year-specific events. It provides a seasonal envelope that supports interpolation and extension. A consequence is that unusual events that deviate strongly from the typical seasonal trajectory can be damped in the reconstructed series, even if they appear clearly in individual clear-sky observations.

Although soil type is related to soil albedo (Post et al., 2000; Rizzo et al., 2023), we excluded it from the RF because of the coarse soil type map (∼ 55.5 km, USDA NRCS, 2023) and due to significant within soil type variation of surface soil color (not shown). Site comparison reveals that the model performs well for soils of intermediate brightness, but its errors increase for very dark or very bright soils (Fig. 3b). This result indicates that the soil type carries information not captured by bulk density, clay, sand, or SOC (e.g., mineral composition and the type of organic matter). It also implies that the importance rankings of predictors influencing α_bs may exhibit spatiotemporal heterogeneity.

Overall, these controls are interacting components governing α_bs rather than independent drivers. Among the intrinsic soil properties, soil texture plays an important role because it establishes the background optical characteristics of the soil surface and regulates soil-water interactions, while soil moisture and solar geometry determine the short-term variability of the observed α_bs.

4.3 Limitations

Despite the process-oriented reconstruction of bare soil albedo across European croplands, several limitations remain, primarily associated with situations in which surface heterogeneity or dynamics operate at spatial or temporal scales finer than those resolved by the available predictors.

First, the mismatched scales between field and satellite observations introduce representativeness errors. α_bs,site is representative for 8–1000 m² (Eichler, 2022) and satellite-based α_bs,RF is for 90 000 m². Soils might show substantial variation in soil properties such as texture, color, crusting, residue, and moisture at the scales of meters (Garten et al., 2007), and the bare soil fraction within a 0.3 km footprint can shift rapidly during sowing and post-harvest. Thus, the resulting α_bs,RF, possibly represents an area-weighted mixture of exposed soil, crop residues, and smooth signals of α_bs dynamics (Yue and Tian, 2020; Dai et al., 2022; Yu et al., 2024).

Second, a sampling bias exists in the training data as periods with crop, clouds or snow covers with extreme soil moisture conditions were excluded. Bare soil is often exposed in autumn or winter periods, only about 35 d in each year from Sentinel-2 data averaged at the 8 ICOS sites satisfied our criteria (Figs. S1 and S2). Fixed local overpass time of Sentinel-2 (∼ 10:30) provides sparse directional sampling of α_bs at the pixel level. Thus, anisotropy is averaged and variance related to extreme cases in the subfield is reduced, with a tendency to underestimate high-magnitude fluctuations (Figs. 2, 5 and 6), such as few observed α_bs,S2 during the periods of 2019 winter at BE-Lon, DE-RuS, and FR-Lam sites (Figs. 6a, b and d) (Sudmanns et al., 2019; Harmel et al., 2018).

Third, base soil masking is not perfect. Crop residue visible or occasional sparse regrowth normally have a higher albedo than the soil (Fig. S5) (Yu et al., 2024), while surface shadows created by terrain or clouds exhibit low reflectance similar to dark soil (Qiu et al., 2019). Using multiple indices in this analysis mitigates but does not eliminate these confusions entirely.

Fourth, the selection of predictors in the RF model might induce uncertainty. The use of coarse resolution (1 km) soil moisture data together with high resolution data could have led to an underestimation of the effect of spatial variation of soil moisture on α_bs. A test of adding soil moisture as a predictor for α_bs,spatial led to the elimination of it during the feature elimination procedure, which deviated from the earlier findings of a well-known negative correlation of spatial patterns of soil moisture and α_bs (Gascoin et al., 2009). Compared with the stable seasonal variation of local solar zenith angle, soil moisture has a higher potential to be disturbed by the environment such as cloud or snow covering periods, which adds noise and reduces the explanatory power. It might explain the higher ranking of local solar zenith angle than soil moisture in predicting α_bs,RF (Fig. 2b and d). In theory, the relative contribution of these two predictors depends on specific scenarios. The comparably higher importance of zenith angle in this study should not be taken to mean that soil moisture is intrinsically less influential in controlling α_bs (Yang et al., 2008; Roosjen et al., 2015). In addition, very sharp post-event shifts of α_bs, especially caused by frequent agricultural activities in sub-regions such as irrigation or tillage, can be muted if the relevant proxy indicators are missing in RF (Mohr et al., 2023; Yu et al., 2024; Goodwell et al., 2025). One important example is varying soil surface roughness which varies according to agriculture practices, rainfall events and crust formation. Its impact on α_bs was not explicitly represented as the required spatiotemporal datasets of soil roughness are currently unavailable at the European scale. In addition, assessing the impact of temporal variation of soil roughness on predicted α_bs at site scale remains challenging as no soil roughness measurements are available at the evaluation sites. Our RF framework relies on intrinsic soil properties that may partially capture roughness-related variability (Herodowicz-Mleczak et al., 2022; Rizzo et al., 2023) but cannot explicitly resolve the effects of agricultural practices, rainfall-induced smoothing, or soil crust formation. In their absence, the model smooths these rapid dynamics and underestimates the abrupt changes in α_bs, yielding a weak response (Fig. 2a and c). Uncertainty of predicted α_bs,RF therefore increases when these rapid physical changes occur on the surface.

Collectively, these scale and sampling properties explain why α_bs,RF preserves seasonal envelopes and regional patterns yet under-represents rapid, high-magnitude extremes. These limitations also highlight the need for future integration of higher-resolution meteorological and soil moisture products, improved surface characterization, and more continuous albedo observations to better capture short-term variability in agricultural landscapes.

5 Conclusion

By integrating explicit bare soil screening of Sentinel-2 surface albedo observation with a random forest reconstruction framework, this study provides a spatially explicit and temporally continuous estimate of bare soil albedo across European croplands at field-relevant scales. The reconstructed dataset captures the pronounced seasonal variability in bare soil albedo and exhibits spatial patterns. Contrary to our initial expectation that soil properties control the spatial distribution, they define relatively stable spatial contrasts in bare soil albedo, without dominating the predictor importance rankings. Additional spatial variability is explained by radiative geometry and temporally averaged soil moisture. The temporal evolution of the estimated soil albedo shows improved coherence with site-scale measurements compared to existing continental satellite products, particularly in resolving seasonal dynamics of soil albedo. This result further indicates that short-term temporal variations are more closely linked to fluctuations in soil moisture than to static soil properties, while variations in solar radiation act as a secondary modulating factor, thereby supporting our second hypothesis regarding the drivers of temporal dynamics.

While uncertainties remain under conditions of rapid surface change such as immediately following rainfall or irrigation events, residue redistribution, or tillage, these limitations primarily reflect constraints in observational sampling and predictor resolution rather than inconsistencies in the reconstructed patterns. Such uncertainties are expected to diminish with the integration of higher-resolution hydrometeorological drivers (e.g., soil moisture) and denser temporal sampling of bare soil conditions. Overall, the dataset is well-suited for characterizing seasonal bare soil albedo variability across heterogeneous agricultural landscapes.

Quantifying high spatiotemporal bare soil albedo provides a physically grounded basis to diagnose the seasonal water and energy dynamics that couple agriculture to climate, especially in agroecosystem modelling. Bare soil albedo directly modulates surface radiation partitioning, influencing soil temperature, evaporation, near-surface energy exchange, and cropland-climate coupling. Incorporation of daily bare soil albedo information into land surface and agroecosystem models has the potential to improve simulations of short-term thermal and moisture dynamics, which are poorly represented when albedo is prescribed at coarse spatial resolution or monthly time steps. Existing modelling approaches typically rely on static or climatological representations of bare soil albedo, limiting their ability to capture rapid transitions that arise from soil wetting-drying cycles and agricultural management. Comparably, the reconstructed albedo presented has more stable temporal evolutions and smoother trends. Incorporating daily bare soil albedo into surface energy balance frameworks enable these dynamics to be explicitly represented, offering new opportunities to diagnose surface energy feedbacks, thereby improving the simulation of temperature stress on plants and soil organisms. Such diagnostics are not attainable with monthly bare soil albedo products due to its temporal aggregation effect. With planned temporal extension and application beyond Europe, this reconstructed bare soil albedo opens new opportunities for climate assessment, crop system evaluation across large spatial and temporal scales, thereby facilitating climate-smart cropland management.