Abstract

Accurate estimation of soil temperature is essential for understanding land–atmosphere interactions and supporting agricultural management in semi-arid regions. This study focuses on the analysis of hourly soil temperature variability and the performance of Multiple Linear Regression (MLR) models at depths of 10, 20, 30, and 40 cm in a tropical semi-arid soil of Burkina Faso. Using in situ meteorological data collected at the Tanghin (Saaba, Kadiogo) experimental site, the results reveal a progressive attenuation of temperature fluctuations with depth, confirming the buffering effect and thermal inertia of subsurface layers. A clear phase lag between soil and air temperatures was observed, increasing with depth and particularly evident beyond 20 cm, which reflects delayed heat propagation within the soil profile. MLR models performed reasonably well at shallow depths, especially when driven by a minimal set of meteorological predictors (Ta, BP, WS, RH), achieving R² values up to 0.75 at 10 cm. However, their predictive ability declined with increasing depth and under wet-season conditions, where soil temperature becomes increasingly decoupled from atmospheric drivers. Moreover, the models struggled to reproduce full diurnal cycles, underestimating nighttime cooling and overestimating during the rainy season, likely due to limited sensitivity to soil moisture dynamics. These findings emphasize both the potential and the limitations of MLR in capturing subsurface thermal variability.

Keywords

Soil Temperature, Subsurface Modeling, Diurnal Variability, Multiple Linear Regression (MLR), Semi-Arid Climate

Share and Cite:

1. Introduction

Soil temperature is a key variable governing the energy and water fluxes at the land-atmosphere interface. It plays a central role in the regulation of evaporation, soil respiration, microbial activity, seed germination, and root development, and thus significantly influences hydrological and agricultural processes [1]-[3]. In dry tropical regions such as Burkina Faso, located in the Sahel, accurate monitoring and forecasting of subsurface soil temperature is essential for sustainable land management, irrigation scheduling, and climate adaptation strategies. Moreover, in the context of increasing climate variability, capturing soil thermal dynamics becomes even more crucial for understanding feedback mechanisms in agro-hydro-meteorological systems [3]. In particular, having access to high-resolution soil temperature data across multiple depths is vital for estimating the components of the surface energy balance, especially the ground heat flux [4] [5]. Furthermore, with the growing frequency of heatwave events in recent years, soil temperature has become a critical variable for studying urban heat island effects and radiative processes influencing near-surface microclimates [6] [7].

Various methods have been used to monitor soil temperature, ranging from in situ sensors to satellite-based remote sensing. Ground-based measurements, typically obtained using thermocouples, thermistors, or advanced sensors, offer high temporal resolution and reliable data near the surface [8] [9]. However, the implementation of such ground monitoring networks remains limited due to high installation and maintenance costs, especially in resource-constrained regions like the Sahel. In many low-income countries, including Burkina Faso, the national meteorological agencies rarely monitor soil temperature as part of their routine observations. Instead, a few agrometeorological institutions attempt to measure it sporadically, often on a limited spatial scale, typically tied to research projects or specific field trials. On the other hand, remote sensing approaches have enabled the retrieval of land surface temperature (LST) and, more recently, subsoil temperature estimates at increasing depths using microwave or thermal infrared sensors. Products such as those from MODIS, SMOS, or SMAP allow soil temperature profiling at regional to global scales [5] [10]. Nevertheless, these satellite-derived datasets often suffer from low spatial resolution and reduced accuracy over heterogeneous terrains, particularly in arid and semi-arid zones like the Sahel. The absence of dense and reliable ground observations further limits the validation and bias correction of satellite products, compromising their usefulness for local-scale applications. This reinforces the need to develop and validate numerical estimation methods that can provide accurate and continuous soil temperature estimates across different depths and environmental conditions. Given these limitations of direct measurement methods, especially in data-scarce regions like the Sahel, it becomes essential to rely on numerical estimation techniques to obtain continuous and spatially representative soil temperature data. In this context, various techniques have been developed to estimate soil temperature.

The most physically grounded approach involves solving the one-dimensional heat conduction equation derived from Fourier’s law, with surface temperature and soil properties as boundary conditions [5] [11]-[13]. While these methods are accurate under well-defined conditions, they require detailed soil thermal properties and often fail to account for complex surface-atmosphere interactions or heterogeneous field conditions [9]. This method is also the simplest among physically based approaches. The more comprehensive and complete involves solving the coupled heat and mass transfer equations, which include both liquid and vapor water flow (using Richards’ equation for unsaturated flow) and heat transfer mechanisms, often under non-equilibrium conditions with phase change between liquid and vapor [14] [15]. Advanced simulation tools such as HYDRUS (2D/3D) implement these coupled processes with high accuracy [16] [17]. However, this modeling framework requires numerous parameters, such as hydraulic conductivity, water retention curves, vapor diffusivity, and thermal conductivity, which are difficult to measure in situ, especially in data-scarce regions where instrumentation is limited and soil characterization is lacking. As a result, the application of such comprehensive models remains challenging in semi-arid and tropical contexts like the Sahel. To overcome these limitations, statistical models such as Multiple Linear Regression (MLR) have traditionally been employed to estimate soil temperature by establishing linear relationships between meteorological predictors and the target variable.

MLR models are computationally efficient and interpretable, making them valuable for applications with limited data or infrastructure. Asadzadeh et al. [8] compared three statistical and machine learning models, SARIMA, Multiple Linear Regression (MLR), and Artificial Neural Networks (ANN), to predict daily soil temperature at multiple depths (5, 10, 20, 50, and 100 cm) at the Riley Station (Oregon, USA). They showed that although MLR performed less accurately than ANN, it provided a simple, robust, and cost-effective modeling approach using easily accessible meteorological variables (e.g., air temperature, infrared radiation). In a similar vein, Cemek et al. [18] assessed various soil temperature estimation techniques (depths: 5 to 100 cm) in the Tokat region, a semi-arid area of Turkey, by comparing MLR with fuzzy logic models (MFIS, ANFIS) and various machine learning algorithms (RF, XGB, MLP, KNN, SVR). MLR was used as a baseline model and, although less accurate than ANFIS or RF, it remained useful for straightforward and interpretable predictions in precision agriculture. In a subtropical grassland during the summer season, Progga et al. [9] evaluated soil temperature at various depths (5, 10, 20, 30, and 50 cm) in a loamy soil under arid tropical conditions in Bangladesh, using meteorological parameters (minimum/maximum air temperatures, relative humidity, pressure, dew point, precipitation, sunshine duration, wind speed, evaporation) through an MLR model. At 5, 10, and 20 cm depths, the MLR model outperformed the HYDRUS-1D model, with R² values of approximately 0.96, 0.94, and 0.88, and MAPE values of 0.81%, 0.87%, and 1.05%, respectively. At 50 cm depth, the R² dropped to 0.55 and MAPE increased to about 1.25%. Their findings suggest that MLR is reliable for modeling soil temperature up to ~20 cm depth. Overall, these studies confirm that MLR remains a robust and interpretable method for estimating soil temperature in subtropical and semi-arid regions, particularly up to ~20 cm depth. More generally, it often serves as a baseline model, being simpler yet less accurate than some advanced algorithms. It was specifically selected for this study due to its interpretability, low computational cost, and straightforward implementation, which are particularly valuable attributes for potential applications in resource-limited contexts.

This study investigates hourly soil temperature variability and estimation at multiple depths (10, 20, 30, and 40 cm) in a tropical semi-arid soil of Burkina Faso using a statistical regression approach (Multiple Linear Regression, MLR) based on in situ meteorological data. While most research emphasizes daily or monthly scales and surface layers, subsurface hourly modeling and diurnal-nocturnal variability remain underexplored. By leveraging high-frequency meteorological observations in a data-scarce tropical context, this work evaluates the performance and limitations of MLR to identify suitable approaches for capturing subsurface thermal dynamics and supporting land surface modeling, precision agriculture, and climate resilience.

2. Materials and Methods

2.1. Site Characteristics and Data Gathering Techniques

The study was conducted in Tanghin, a village located in the rural municipality of Saaba in the Kadiogo Province of Burkina Faso (12˚22’5"N, 1˚25’0"W). The site lies in the Sudano-Sahelian climatic zone, characterized by two distinct seasons: a dry season from November to May and a rainy season from June to October. The average annual rainfall ranges from 600 to 900 mm, and the average temperature varies between 22.6˚C and 35.5˚C [19].

For ground data collection, an excavation was made to a depth of 80 cm to facilitate sensor installation and ensure accurate soil measurements. Within this soil profile, two distinct layers were identified, a subsurface soil up de 30 cm depth. The physical properties of the soil at 10 - 30 cm and 30 - 80 cm depths are summarized in Table 1. Both layers exhibit a sandy loam texture, with nearly identical porosity (≈0.41) and bulk density (~1560 kg·m⁻³), suggesting homogeneity in the soil’s water retention and thermal characteristics within the measured horizon. Vegetation was periodically cleared to minimize interference with surface sensors, particularly during the rainy season.

Table 1. Physical soil properties by profile.

The meteorological and soil data measurements were collected using instruments from Campbell Scientific Inc. (Logan, UT, USA) and METER Group. The ClimaVUE50 autonomous weather station was used to record air temperature, relative humidity, vapour pressure, barometric pressure, wind speed and direction, solar radiation, sunshine duration, precipitation, and lightning activity. A CR350-CELL215 datalogger facilitated data collection and transmission via a 4G modem. Additional atmospheric parameters derived from direct measurements include Dew Point Temperature, Wet Bulb Temperature, Actual Vapor Pressure, Saturation Vapor Pressure, Vapour Pressure Deficit (VPD), and Evapotranspiration were calculated from raw measurements using internal algorithms. To monitor some soil data, TEROS21 and CS655 sensors were installed at multiple depths (10, 20, 30, and 40 cm). These instruments enabled the continuous collection of soil temperature, volumetric water content, and matric potential. All variables were logged at three temporal resolutions: every 10 minutes, hourly, and daily. Figure 1 displays the station setup.

Figure 1. Field instrumentation used for soil and meteorological data acquisition at the Tanghin experimental site (Saaba, Burkina Faso). (Left) The ClimaVUE50 autonomous weather station, mounted on a metallic frame and powered by a solar panel system, records key atmospheric variables. (Center) Internal view of the control unit containing the CR350-CELL215 datalogger and wiring for real-time data transmission. (Right) Cross-sectional schematic of the sensor layout showing TEROS21 & CS655 probes deployed at 10, 20, 30, and 40 cm depths, each individually connected to the datalogger.

For the present study, we selected hourly data for the following meteorological variables, denoted by their abbreviations: air temperature (Ta, ˚C), relative humidity (RH, %), dew point temperature (DPT, ˚C), wind speed (WS, m/s), wind direction (WD, ˚), barometric pressure (BP, hPa), precipitation (RAIN, mm), solar radiation (RAD, W/m²), and vapour pressure deficit (VPD, hPa). Soil temperature at various depths is denoted by STxy, where x indicates the depth in cm (10, 20, 30, 40) and y specifies the sensor type: t for TEROS21 and c for CS655 (e.g., ST10t = soil temperature at 10 cm from TEROS21).

The data used in this study were collected from 22 September 2023 to 08 March 2025, covering a total period of 534 days. Measurement instrument failures occurred between 01 April 2024 and 14 May 2024, and between 03 and 28 July 2024, resulting in missing data. In total, the dataset comprises 464 days and 11,136 hours (corresponding to 11,136 data rows) for this study. Missing data entries were deleted prior to analysis. Indeed, since the MLR model itself does not account for or rely on any temporal structure, this deletion method cannot introduce a time-dependent bias in the model’s estimates or performance.

2.2. Methodology

To estimate soil temperature from in situ meteorological inputs, we adopted a statistical approach based on multivariate linear regression (MLR). This classical method is widely used to model the linear relationship between a dependent variable, soil temperature in this study, and several independent predictors [8]. MLR relies on the least squares principle to minimize the sum of squared errors between observed and predicted values. For the model to yield reliable and interpretable coefficients, it is crucial that the predictors exhibit low intercorrelation, thereby avoiding redundancy and reducing the variance of estimators.

2.2.1. Variable Selection and Multicollinearity in MLR Baseline Model

To handle potential multicollinearity among the predictors [20], a two-step procedure was implemented. First, a stepwise variable selection based on the Akaike Information Criterion (stepAIC) was used to retain the most relevant variables. Then, the Variance Inflation Factor (VIF) was computed to quantify the extent of multicollinearity. Variables with VIF > 5 were treated with caution, while those exceeding a threshold of 10 were considered highly collinear; in such cases, only one representative variable was retained [9] [21]. This strategy ensures a parsimonious and stable regression model suitable for explanatory and predictive purposes.

Figure 2. Variance Inflation Factor (VIF) values for soil temperature at different depths with respect to the meteorological predictors.

Figure 2 displays the VIF values associated with each meteorological variable included in the multivariate linear regression models for soil temperature at four depths (ST10t to ST40t). The variables VPD (Vapour Pressure Deficit) and RH (Relative Humidity) exhibit very high VIFs (>10), indicating strong multicollinearity, especially for VPD (~17.9), which could severely bias regression estimates. Air temperature (Ta) also shows moderate multicollinearity (VIF ≈ 6.36). To maintain model reliability, variables with high VIFs were either excluded or carefully selected based on correlation structure and predictive relevance.

2.2.2. Development of Regression Scenario

To apply multivariate linear regression (MLR) across all soil depths, we established a set of input scenarios based on combinations of meteorological predictors, following a strategy similar to that proposed [3] for tree-based hybrid Data Mining models. For each depth, air temperature (Ta), a variable that is easy to measure and broadly available, is used as the baseline predictor. Additional variables were progressively added to account for measurement complexity and to minimize multicollinearity, assessed using Variance Inflation Factor (VIF) values after an initial selection based on stepwise AIC (stepAIC). For example, in resource-limited settings with fewer ground-based instruments, relative humidity (RH) may be favored over vapor pressure deficit (VPD), which requires two separate measurements.

Based on this approach, six (06) regression scenarios were designed to support the implementation of MLR models. These scenarios reflect a balance between data availability, predictor relevance, and acceptable levels of multicollinearity. Identical sets of predictors were retained for the shallow soil layers (10 and 20 cm) and, similarly, for the deeper layers (30 and 40 cm), following the outcomes of stepAIC selection and VIF screening. Table 2 and Table 3 summarize these scenarios and provide the rationale for each variable’s inclusion. The two tables differ in the final variable included: S6 adds rainfall (RAIN) for shallow layers (ST10 - ST20), while it adds wind direction (WD) for deeper layers (ST30 - ST40), reflecting the differing relevance of surface vs. subsurface meteorological influences.

Table 2. Description and rationale of regression scenarios applied to shallow soil layers (10 cm and 20 cm depths).

Table 3. Description and rationale of regression scenarios applied to deeper soil layers (30 cm and 40 cm depths).

2.3. Criteria for Model Evaluation

To evaluate and compare the performance of the different MLR scenarios tested in this study, we used a set of classical regression metrics, including the coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE). Additional information criteria such as the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) were computed to assess model parsimony and penalize overfitting based on the number of input variables retained. These indicators quantify both the explanatory power and the predictive accuracy of the models across soil depths

The formulation of each evaluation metric is presented as follows:

• Mean Absolute Error (MAE)

$\text{MAE=}\frac{{\displaystyle {\sum }_{i=1}^N\left|y_i-{\widehat{y}}_i\right|}}{N}$ (1)

• Coefficient of Determination (R²)

${\text{R}}^2=1-\frac{{\displaystyle {\sum }_{i=1}^N{\left(y_i-{\widehat{y}}_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^N{\left(y_i-\overline{y}\right)}^2}}$ (2)

• Root Mean Square Error (RMSE)

$\text{RMSE=}\sqrt{\frac{{\displaystyle {\sum }_{n=1}^N{\left(y_i-{\widehat{y}}_i\right)}^2}}{N}}$ (3)

where N is the total number of observations, y_i is the observed value at time i, ${\widehat{y}}_i$ is the predicted value, and $\overline{y}$ is the mean of the observed values.

• Akaike Information Criterion (AIC)

$\text{AIC}=-2\cdot \log \left(L\right)+2k$ (4)

• Bayesian Information Criterion (BIC)

$\text{BIC}=-2\cdot \log \left(L\right)+\log \left(n\right)\cdot k$ (5)

where n is the number of observations, k is the number of model parameters, and L is the maximized value of the likelihood function of the model.

3. Results and Discussion

3.1. Overview of Soil Temperature Variability

Figures 3-5 present the variability of hourly soil temperature measurements during the study period. Figure 3 shows the overall temporal evolution of the different depths soil temperature measured by TEROS21 across the entire period. It highlights hourly fluctuations and missing data intervals (shaded areas), as an illustration among other recorded values. A higher level of noise is observed near the surface, while measurements beyond 20 cm appear smoother. This is due to instantaneous variability in surface-layer processes such as incoming shortwave radiation, surface-emitted longwave radiation, and rapid fluctuations in atmospheric conditions (wind, air temperature, humidity).

Figure 3. Hourly Soil Temperatures at the four depths with daily moving average and missing data intervals (shaded areas).

Figure 4. Multi-day diurnal dynamics of soil and air temperatures (15-18 June 2024).

Figure 5. Monthly Mean Soil Temperatures in 2024 by depth.

Figure 4 reveals a clear attenuation of temperature fluctuations with increasing soil depth, highlighting the buffering role of deeper layers against diurnal variations. These results show that the soil temperature is higher than the air temperature at all depths, although this difference is more pronounced in the near surface layers. The comparison with air temperature confirms the soil’s thermal inertia, particularly beyond 20 cm. Interestingly, during the morning hours, soil temperature continues to decrease while air temperature begins to rise. For the 10 cm depth, the timing of minimum temperature lags behind air temperature by more than two hours, as annotated in Figure 4. This lag increases with depth; for instance, at 40 cm, the soil temperature reaches its minimum almost simultaneously with the peak in air temperature, highlighting the strong thermal inertia of deeper layers during the rainy season.

This shift likely results from the delay in heat propagation through the soil profile, governed by radiative fluxes and energy exchanges at the surface. A more detailed investigation is needed to understand how this phase shift varies across the rainy and dry seasons, particularly in the context of arid tropical environments. Conversely, maximum temperatures occur almost simultaneously across layers, suggesting faster surface heating during daytime.

Figure 5 summarizes the monthly mean soil temperatures using a heatmap, providing insight into seasonal variations with depth.

These visualizations provide insight into vertical and seasonal thermal fluctuations and introduce the subsequent analyses focused on distribution patterns using violin plots and correlations with explanatory variables.

3.2. Exploration of Meteorological and Soil Temperature Distributions

To further explore the distribution and variability of the key environmental variables recorded at the Saaba station, violin plots were generated for air and soil temperatures as well as major meteorological parameters.

Figure 6. Distribution of meteorological and soil variables over the observation period. The top panel shows air temperature (Ta), dew point temperature (DPT), and soil temperature (ST) at 10, 20, 30, and 40 cm depths measured with TEROS21 sensors (t) and CS655 sensors (c). The bottom panels illustrate (i) wind speed (WS) and wind direction (WD), (ii) relative humidity (RH) and barometric pressure (BP), and (iii) vapor pressure deficit (VPD) and solar radiation (RAD). Mean values are displayed for each distribution.

Figure 6 provides visual summaries of the central tendencies and variability of the environmental variables over the observation period. The top panel shows the distributions of air temperature (Ta), dew point temperature (DPT), and soil temperatures (ST) at 10, 20, 30, and 40 cm depths, measured with TEROS21 (t) and CS655 (c) sensors. Air and soil temperatures exhibit more stable distributions. Air temperature (Ta) averages 26.8˚C, while dew point temperature (DPT) is significantly lower at 13.3˚C, indicating a marked vapor pressure gap. Soil temperatures are more stable and consistently average around 31˚C - 32˚C across depths and sensor types, with slight increases in the TEROS21 sensors at 10 and 30 cm. The close agreement between TEROS21 and CS655 measurements supports their interchangeability, although TEROS21 data will be preferred for further modeling.

The bottom panels display the distributions of key meteorological drivers. Wind speed (WS) shows a narrow distribution centered around 0.8 m/s, suggesting weak air movement, whereas wind direction (WD) is broadly distributed with a central tendency around 151.9˚, reflecting variable wind patterns. Relative humidity (RH) exhibits a bimodal distribution with substantial variability (mean ≈ 51.5%), while barometric pressure (BP) follows a more symmetrical distribution centered near 978.8 hPa. Vapor pressure deficit (VPD) spans a wide range with a positively skewed distribution, indicating frequent episodes of atmospheric dryness. Solar radiation (RAD) also shows a right-skewed pattern, with a mean of approximately 217.6 W/m², reflecting the expected diurnal variability in tropical settings.

3.3. Correlation Analyses between Meteorological Predictors and Soil Temperatures

To complement the distribution analysis, a correlation matrix was constructed to explore the relationships between meteorological variables and soil temperatures at multiple depths using the two types of sensors (TEROS21 and CS655). Figure 7 summarizes the pairwise Pearson correlation coefficients, offering insights into the interactions driving subsurface thermal dynamics.

The correlation matrix reveals strong positive associations between air temperature (Ta) and soil temperatures across all depths, especially at shallower levels (e.g., ST10t: r = 0.83), highlighting the dominant influence of surface conditions. Other meteorological variables such as vapor pressure deficit (VPD) and solar radiation (RAD) also show moderate to weak positive correlations, suggesting their role in surface energy balance.

Barometric pressure (BP), in contrast, exhibits strong negative correlations with soil temperature, indicating an inverse relationship possibly linked to broader atmospheric processes rather than local heating. Dew point temperature (DPT) gains relevance with increasing depth, surpassing Ta in explanatory strength beyond 30 cm.

The soil temperature variables demonstrate remarkably high internal consistency, both within and between sensor types (e.g., r > 0.98 between ST30c and ST40c, and r = 0.99 between ST10t and ST10c), confirming the robustness of measurements and justifying the use of TEROS21 data for subsequent analyses. These results underline the stability of vertical soil temperature profiles and support the representativeness of selected metrics.

Figure 7. Pearson correlation matrix between meteorological variables and soil temperatures measured by TEROS21 (ST10t to ST40t) and CS655 (ST10c to ST40c) sensors.

To further dissect these relationships, Figure 8 presents the Pearson correlation coefficients between soil temperature at four depths (ST10t to ST40t, measured by TEROS21) and key meteorological variables. This bar plot format allows for a clearer view of how the strength and direction of associations evolve with soil depth. As observed, air temperature (Ta) and barometric pressure (BP) remain the most influential factors across all depths. Ta shows a strong positive correlation that progressively weakens with depth (from r = 0.83 at 10 cm to r = 0.58 at 40 cm), highlighting the decreasing influence of surface heating. Conversely, BP maintains a consistently strong negative correlation (r < –0.65), suggesting broader atmospheric control mechanisms. A plausible physical explanation is that high-pressure systems are typically associated with clearer skies, which enhance nocturnal radiative cooling of the land surface, and with the advection of cooler, drier air masses. This combination of factors leads to a more pronounced drop in soil temperature, hence the observed strong negative correlation. Dew point temperature (DPT) gains importance at deeper layers, surpassing Ta from 30 cm downward, thus reflecting latent heat and moisture effects in subsurface processes. Other meteorological variables, including vapor pressure deficit (VPD), relative humidity (RH), and wind parameters, exhibit weak to moderate correlations, with no clear patterns across depths. These findings emphasize the progressive dampening of meteorological influence with increasing soil depth, and reinforce the selection of key predictors in subsequent modeling efforts.

Figure 8. Pearson correlation coefficients between soil temperature at 10, 20, 30, and 40 cm (measured by TEROS21 sensors) and selected meteorological variables.

3.4. Comparison of Observed and MLR-Predicted Distributions of Hourly Soil Temperature Across Depths

The results obtained for each regression scenario are presented using violin plots, which allow simultaneous comparison of both the distribution shape and the central tendency (mean and median) of the measured and estimated soil temperatures. This visual representation helps assess the accuracy and consistency of the multivariate linear regression models across different soil depths. Then, Figure 9 compare the observed and predicted soil temperatures across all the four depths using the six MLR scenarios.

The distribution shapes offer more insight. At shallow depths (10 and 20 cm), adding variables slightly improves the resemblance between predicted and observed distributions, reducing the asymmetry and narrowing the spread. At deeper layers (30 and 40 cm), the gain is less evident, and the predicted distributions remain narrower than the measured ones, suggesting that the models capture the central tendency well but may underestimate the full variability. Overall, increasing the number of predictors improves the representation of the distribution’s shape, particularly at surface levels, but has diminishing returns at depth due to the smoothing effect of the soil profile.

3.5. Performance Evaluation of Tested MLR Models Using Scatter Plots, Regression Metrics, and Fitted Equations

To evaluate the predictive performance of the multivariate linear regression (MLR) scenarios, scatterplots of predicted versus observed soil temperatures were generated for each scenario and soil depth. Key performance metrics, including the coefficient of determination (R²), root mean square error (RMSE), mean absolute error (MAE), and model equations, were used to assess model accuracy and goodness of fit.

Figure 9. Distribution comparison of observed vs. predicted hourly soil temperature using MLR models under six scenarios, at four soil depths: (a) 10 cm; (b) 20 cm (c) 30 cm and (d) 40 cm.

Figures 10-13 summarize the results across the six scenarios for estimating hourly soil temperature at four depths (10, 20, 30, and 40 cm), by comparing observed and predicted values in each case by scatterplot.

For 10 cm-depth (ST10), performance steadily improves with model complexity. The R² increases from 0.645 (S1) to 0.745 (S6), while RMSE decreases from 2.77˚C to 2.34˚C. This indicates that adding predictors like wind speed, pressure, and solar radiation enhances the fit without overfitting. The S6 scenario achieves the best compromise between accuracy and model complexity. At 20 cm-depth (ST20), the models also improve as predictors are added. R² rises from 0.301 (S1) to 0.627 (S6), and RMSE drops from 2.81˚C to 2.05˚C. Notably, adding radiation and wind direction from scenario S5 onwards significantly boosts performance, suggesting greater influence of these factors at this depth.

For 30 cm-depth (ST30), predictive power remains moderate. R² starts low at 0.109 (S1) but reaches 0.475 (S6). RMSE decreases from 2.60˚C to 2.01˚C. Although the accuracy improves with more variables, the dispersion remains more compact, indicating that deeper soil temperature dynamics are harder to capture using surface meteorological inputs alone. At 40 cm-depth (ST40), similar trends are observed. The basic model (S1) yields a low R² of 0.077 and RMSE of 2.38˚C, while the most complex scenario (S6) improves to R² = 0.361 and RMSE = 1.98˚C. However, gains from S4 to S6 are marginal, highlighting the diminishing influence of surface variables with increasing depth.

Figure 10. Observed vs. predicted hourly soil temperature using MLR models under six scenarios at 10 cm soil depth.

Figure 11. Observed vs. predicted hourly soil temperature using MLR models under six scenarios at 20 cm soil depth.

Figure 12. Observed vs. predicted hourly soil temperature using MLR models under six scenarios at 30 cm soil depth.

Figure 13. Observed vs. predicted hourly soil temperature using MLR models under six scenarios at 40 cm soil depth.

Figure 14. Performance metrics (R², RMSE, MAE) of multivariate linear regression models across six scenarios and four soil depths (10, 20, 30, and 40 cm).

Figure 14 presents the performance metrics, coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE), for the MLR models across the six scenarios and the four soil depths. These results support the scatterplot analysis, highlighting that the inclusion of additional meteorological predictors improves model accuracy, especially for the surface layers (10 - 20 cm), where R² exceeds 0.6 in the most complete scenarios while maintaining acceptable error levels (RMSE < 2˚C, MAE < 1.7˚C).

We summarize in Tables 4-7 the performance metrics of the MLR models across six scenarios for estimating hourly soil temperatures at four different depths (10, 20, 30, and 40 cm). For each scenario, we add also Akaike information criterion (AIC), Bayesian information criterion (BIC), and the fitted regression equation.

Table 4. Performance metrics and model equations for MLR across six scenarios for predicting hourly soil temperature at 10 cm depth.

Table 5. Performance metrics and model equations for MLR across six scenarios for predicting hourly soil temperature at 20 cm depth.

Table 6. Performance metrics and model equations for MLR across six scenarios for predicting hourly soil temperature at 30 cm depth.

Table 7. Performance metrics and model equations for MLR across six scenarios for predicting hourly soil temperature at 40 cm depth.

Across all depths, the AIC and BIC values decrease gradually from Scenario S1 to S6, indicating a better model fit as more explanatory variables are added. However, the variations in AIC and BIC between successive scenarios, especially between S4, S5, and S6, are relatively modest. This implies that while model accuracy slightly improves, the added complexity may not always be justified, particularly for operational use in data-sparse regions. Intermediate models (such as S3 or S4) may offer an optimal trade-off between performance and parsimony.

Overall, and to summarize the MLR models in estimating hourly soil temperatures across depths, we can say that predictive accuracy is strongly influenced by the number and relevance of meteorological predictors. At shallow depths (10 cm and 20 cm), the models show substantial improvement as more variables are included, with R² values increasing from ~0.65 (baseline with Ta alone) to ~0.75 (for ST10) and ~0.63 (for ST20) in the most complete scenarios. In contrast, deeper layers (30 cm and 40 cm) remain less predictable with R² values not exceeding 0.48 and 0.36, respectively, regardless of model complexity. This confirms that subsurface temperatures are more insulated from short-term atmospheric variability.

Our results are consistent with those of [9] who also found that at soil depths greater than 30 cm, both models failed to accurately estimate soil temperature using daily data. Either method, however, could satisfactorily predict soil temperature at depths between 0 and 20 cm. In their study, model performances were notably higher than in ours, with R² values of 0.96, 0.94, and 0.88 at 5, 10, and 20 cm depths, respectively, using daily data and 13 input variables including minima and maxima. By contrast, our hourly model yielded R² values of 0.745 and 0.627 at 10 and 20 cm, respectively, using only seven meteorological parameters with multicollinearity excluded. This difference in performance may be attributed to the higher noise level in hourly data compared to daily data.

So, for practical estimation with regression methods of soil skin temperature (ST10) at hourly scale with a good approximation (targeting R² ≈ 0.75), a minimum set of four variables is sufficient: air temperature (Ta), barometric pressure (BP), wind speed (WS) and relative humidity (RH). These variables, commonly measured in standard weather stations, capture the essential surface energy exchanges that drive soil temperature variations. Adding further variables like solar radiation (RAD) and wind direction (WD) slightly improves model precision.

3.6. Temporal Consistency and Variability Assessment of Predictions in Soil Temperature Series

While standard performance metrics (e.g., R², RMSE, MAE), scatter plots, and distribution-based visualizations such as violin plots offer valuable insights into model accuracy and global behavior, they often fail to reveal local temporal inconsistencies in time series data. In the case of soil temperature estimation, which involves spatio-temporal sequences governed by cyclic thermal processes, it is crucial to verify whether the predicted values follow the natural variability and dynamics of the observed series, including the correct reproduction of local extrema (minima and maxima) across different depths. This aspect is particularly important from a physical standpoint, as soil temperature responds to radiative forcing, surface energy balance, and subsurface heat propagation, which collectively induce characteristic oscillations in the signal. A model that performs well globally but fails to replicate these key temporal patterns may lead to misleading conclusions in process-based applications. Therefore, we systematically compare the temporal structure of the predicted and observed series with a specific focus on the preservation of diurnal and sub-daily variability. This temporal fidelity analysis complements the global evaluation metrics and ensures that models not only approximate values accurately but also respect the physical dynamics underlying soil temperature variations.

So, we present in Figure 15 the temporal variability of measured and predicted soil temperatures under different MLR scenarios (S1, S4, S6), along with their 24-hour rolling means and 95% confidence intervals. Scenario S1 was selected because it relies solely on air temperature (Ta), a variable that is easily and widely measured. Scenario S4 incorporates additional commonly available meteorological variables: wind speed (WS), barometric pressure (BP), and relative humidity (RH), while scenario S6 includes all meteorological predictors, with variable selection adjusted to account for multicollinearity. The analysis focuses exclusively on the 10 cm soil depth, where the MLR predictions yielded the highest statistical performance (R² > 0.60 across all scenarios), unlike deeper layers where the predictive accuracy of MLR models was markedly lower.

The MLR-based predictions demonstrate good agreement with observations during the dry season (January-February-March, JFM), but tend to overestimate soil temperature during the rainy season (June-September) and underestimate it toward the end of the wet season (October-November-December, OND). These deviations suggest that the models lack sensitivity to soil moisture dynamics, a key variable influencing thermal behavior through evaporative cooling and heat capacity modulation. Specifically, higher moisture content increases the soil’s thermal conductivity, promoting heat transfer, while simultaneously raising its heat capacity, requiring more energy to change its temperature. The absence of this variable likely explains the model’s error in simulating the amplitude and phase of diurnal temperature cycles, particularly during wet periods. Incorporating soil water content as an additional predictor could therefore enhance model performance, especially under conditions of high atmospheric humidity and rainfall variability, as highlighted in Taheri et al. [10].

Figure 15. Temporal variability of observed and predicted soil temperature at 10 cm depth under three MLR scenarios (S1, S4, S6), with 24-hour rolling averages and 95% confidence intervals.

To further investigate the ability of the models to reproduce fine-scale temporal fluctuations, we conducted a focused comparison over a short time window (15-18 June 2024), which covers a full diurnal cycle over three days. Figure 16 presents the observed and predicted soil temperatures at 10 cm depth under three MLR scenarios (S1, S4, S6), during 15-18 June 2024, with emphasis on the timing and amplitude of daily extrema.

The short-term validation highlights clear performance differences between the modeling approaches. The predictions systematically underestimate the observed temperatures, especially during the rainy season (June), and fail to capture the full diurnal cycle. This further confirms the limited ability of MLR models to generalize under high humidity and soil moisture conditions typical of the wet season.

Figure 16. Short-term validation of multiple linear regression (MLR) predictions at 10 cm depth for three scenarios (S1, S4, S6) from 15 to 18 June 2024.

4. Conclusions

This study presents a comprehensive comparative analysis of statistical MLR approaches to estimate hourly soil temperature at different depths of 10, 20, 30, and 40 cm in a semi-arid region of Burkina Faso using in situ meteorological data. Conducted at the Tanghin (Saaba, Kadiogo) experimental site, the study highlights the challenges and opportunities associated with modeling subsurface thermal dynamics in low-instrumentation contexts and under high temporal resolution.

Analysis of the variability further revealed a clear attenuation of temperature fluctuations with depth, confirming the buffering effect and thermal inertia of subsurface layers. The observed phase lag between soil and air temperatures, which increases with depth, particularly beyond 20 cm, emphasizes the delayed heat propagation within the soil profile. A more detailed investigation is needed to understand how this phase shift varies across the rainy and dry seasons, particularly in the context of arid tropical environments.

MLR models demonstrate acceptable performance at shallow depths when fed with a minimal set of meteorological variables (Ta, BP, WS, RH), achieving R² values up to 0.75 at 10 cm. However, their predictive ability sharply declines in deeper layers, where the thermal signal becomes increasingly decoupled from surface atmospheric forcing. In addition, the models struggled to accurately reproduce the full diurnal cycles, as evidenced by a consistent underestimation of the diurnal temperature range. This is manifested through an underestimation of nighttime cooling and, to a lesser extent, daytime warming. Furthermore, the models systematically underestimated temperatures during the rainy season. These combined biases reflect a limited sensitivity to soil moisture dynamics, which critically modulate thermal behavior through evaporative cooling and heat capacity effects.

Despite its scope, this study has limitations. Data from a single site may restrict transferability to other contexts with different soil textures, vegetation cover, or climatic regimes. Furthermore, the absence of continuous soil moisture measurements limited the ability to account for its critical role in subsurface thermal regulation.

Future research should focus on extending these models across multiple agroecological zones, integrating soil moisture and thermal properties into prediction frameworks, and exploring hybrid physics-informed machine learning methods to enhance robustness and interpretability. Coupling soil temperature modeling with crop growth, irrigation scheduling, and land surface models would provide valuable tools for precision agriculture and climate resilience in semi-arid regions.

Acknowledgements

The authors have reviewed and edited the output and take full responsibility for the content of this publication.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References