Estimating high-resolution soil moisture by combining data from a sparse network of soil moisture sensors and remotely sensed MODIS LST information

Soil moisture (SM) is a crucial hydrological variable for many aspects of environmental processes and human life. Consequently, it is an indispensable variable for crop growth as it describes the water availability for plants and also constitutes a pivotal variable as it controls the energy exchanges between the soil and the atmosphere, as well as the water movement from the surface to the subsurface (Fang & Lakshmi 2014). Improved knowledge of the SM regime and its interrelationship with the various hydrological processes will facilitate the improvement of meteorological, climatological and hydrological forecasts and the more realistic allocation of water resources (Pablos et al. 2016a; Kofidou & Gemitzi 2023). Furthermore, the enhanced characterization of SM variability at local scales is of paramount importance for hydrological and climate studies (Vergopolan et al. 2022).

A considerable number of studies have been carried out with the objective of acquiring SM at high spatial resolution. These studies have primarily focused on the downscaling of SM from passive microwave radiometer sensors on board satellite missions such as Soil Moisture Active Passive (SMAP) and Soil Moisture Ocean Salinity (SMOS), at the global or local scales. In those works, ground measurements, when available, are employed to validate the results of the downscaling from various downscaling techniques or to bias correct the remotely sensed SM (Kornelsen et al. 2015). The majority of methodologies employed to downscale satellite-derived SM rely on the thermal inertia relationship between surface temperature and soil water content. Such a coupling is well known and has been documented in numerous previous works (Minacapilli et al. 2009; Pablos et al. 2016a; Fang et al. 2018a, 2021a) indicating also that it is strongly dependent on the accurate estimation of the soil thermal conductivity (Minacapilli et al. 2009). The surface brightness temperature is also influenced by the physical temperature of the soil and the surface water content and the vegetation conditions (Pablos et al. 2016a), and consequently, the satellite missions that monitor LST along with other key environmental variables, providing Thermal Infrared (TIR) remotely sensed data, have been extensively employed to acquire SM at the fine resolution (Hain et al. 2011; Piles et al. 2014; Pablos et al. 2016a; Fang et al. 2018a, 2018b, 2021a).

Although the downscaling of satellite-derived SM can be beneficial in areas with no in situ SM monitoring (Fang et al. 2021b), previous works have also highlighted the potential to produce high-resolution SM from a sparse network of SM sensors (Ochsner et al. 2019). The computation of gridded SM from a sparse network of SM sensors presents a significant challenge. The typical spacing of these sensors does not permit the direct application of interpolation techniques, using ground SM observations from sensors. Thus, the possibility of obtaining high-resolution SM datasets, for example, with a resolution of approximately 1 km, is limited. Therefore, other variables are usually evaluated for their predictive ability in relation to SM through the application of a regression formula or other modeling approaches. To illustrate, the Oklahoma Mesonet environmental monitoring network described in Ochsner et al. (2019) encompasses 120 monitoring points with an average distance of 30 km between them. In order to acquire gridded daily SM at 800 m spatial resolution from the Oklahoma Mesonet multi-depth monitoring points, a multi-linear regression method using as predictor variables sand content and antecedent precipitation index (API), was developed (Ochsner et al. 2019). The sand content was available from the Gridded Soil Survey Geographic Database (Soil Survey Staff 2014) at 10-m resolution, whereas the API was estimated from precipitation estimates at 4-km resolution. In a subsequent step, residuals of the measured and estimated SM values were evaluated and an empirical semivariogram of the residuals was computed. Finally, daily SM maps were generated for each examined depth, i.e., 5, 25, and 60 cm, summing the regression model values with the interpolated residuals. Therein, the authors concluded that their approach is suitable to describe the mesoscale SM patterns using a sparse, large-scale ground monitoring network. Nevertheless, there is still a gap in the available techniques when there is no large-scale monitoring network and a limited number of SM sensors operate within a geographically constrained area, for instance, at the provincial level, particularly when there is a lack of detailed information regarding soil types.

In the present work, we exploit the SM and LST relationship to generate gridded daily SM at fine spatial resolution, i.e., 1 km using a sparse network of SM sensors. We argue that different locations (altitude, soil type, land cover characteristics) have different relationships with LST, so we extracted location-based SM–LST relationships and their coefficients were interpolated over the study area, thus providing a novel method for obtaining high-resolution SM. Our approach has certain advantages, as it does not require a large monitoring network of SM and exploits the global availability of remotely sensed LST, merging in situ information with satellite-based datasets.

The dataset consists of daily averages of SM in six in situ monitoring points operated from 1 October 2021. The three monitoring points in Kerasia, Komotini and Nimfea are operated by Democritus University of Thrace (DUTH), equipped with Sentek Drill and Drop single-point SM and soil temperature sensors (https://sentektechnologies.com/products/soil-data-probes/drill-drop/, accessed 28/11/2023), which are installed at 5 cm depth, providing SM measurements in volumetric water content (%), i.e., [volume of water (m³)/volume of soil (m³)] × 100. Those three sensors were used for model development. The two sensors, i.e., Asomati and Karidia, are provided by local farmers operating Watermark SM sensors at 5 cm depth as well, as part of the Davis Pro2 meteorological stations (https://www.davisinstruments.com/products/soil-moisture-sensor-vantage-pro-and-vantage-pro2, accessed 28/11/2023), whereas the same SM sensor has been installed on 5, August 2023 in Ifestos village as part of the DUTH's SM monitoring program, providing SM data until today. Those three sensors were used for model testing. In Table 1, the operation period and the basic annual statistics of the examined variables, along with altitude and soil types in the monitoring stations in the Komotini test area, are provided.

To unveil the possible relation of daily SM with LST, an analysis was conducted of the correlation of in situ SM from the monitoring network and remotely sensed LST from TIR sensors. There are many satellite missions that monitor LST. LST products from Terra and Aqua satellites carrying the MODIS sensor provide daytime and nighttime LST at 1-km spatial resolution at the daily timestep. The algorithms to derive LST products are comprehensively evaluated and are validated in a wide range of environmental conditions (Fang & Lakshmi 2014; Eleftheriou et al. 2018; Lu et al. 2018). Therefore, the latest collection, i.e., collection 6.1 of the daily daytime/nighttime MODIS LST products (MYD11A1v061 and MOD11A1v061) (Wan et al. 2021), was evaluated for possible predictive ability for daily averaged SM. The correlation analysis conducted comprised MODIS Terra and Aqua daytime and nighttime LST, and the associated diurnal temperature range, i.e., the difference between the daytime and nighttime LST of Terra and Aqua satellites, respectively.

The theoretical background of the thermal inertia-based method for the estimation of SM can be found in (Cheruy et al. 2017; Matsushima 2019). Thermal inertia, which describes the resistance of soil to temperature changes, is a property of soil and it is defined as the square root of the product of the volumetric heat capacity of soil and the thermal conductivity, both increasing with the increasing SM. Thermal inertia is found in the formulation of soil heat flux when LST is the only known parameter, leading to the description of land surface processes using satellite-derived LST, and it has been shown that thermal inertia and consequently LST can be used in the estimation of SM (Matsushima 2019).

Equation (1) estimates daily SM (%) at the monitoring station i and at the t date, where is the slope of the regression equation and is the intercept at each monitoring point. is the daytime Aqua LST (^oC) at 1-km spatial resolution at daily timestep, for the specific pixels where monitoring stations are located.

In Equation (3) (%) stands for the corrected daily SM at 1 km, (%) are the interpolated daily residuals at 1 km and (%) is the daily SM at 1 km from Equation (1). Please note that i represents in situ monitoring points in Equations (1) and (2), whereas it stands for individual 1-km grid cells in Equation (3).

The acquired residual corrected daily SM estimates at 1 km were then verified from the three verification in situ monitoring points in Asomati, Kerasia and Ifestos, computing the performance indicators stated below:

In addition, to evaluate the performance of the methodology presented herein with other methods of acquiring high-resolution daily SM, a comparison against (a) simple interpolation and (b) downscaled 1-km daily SMAP (Fang et al. 2022), is presented and discussed in the Results and Discussion sections. The downscaled SMAP SM was developed based on regression equations estimated between SM and LST difference from the Global Land Data Assimilation System–Noah Land Surface Model in 10 NDVI groups from the Advanced Very High-Resolution Radiometer (AVHRR).

Correlation analysis of in situ SM with the Aqua and Terra daytime and nighttime LST provided the results shown in Table 2, indicating that Aqua daytime LST is the parameter with the highest anticorrelation to daily averaged SM in all six monitoring stations.

Analogous findings were indicated in previous works, attributing them to the better capturing of the topographic effects on LST in the daytime MODIS passes, whereas Aqua daytime passes (ascending pass with equator crossing time at 13:30) are closer to the daily maximum LST, which is considered a better variable to downscale SM (Pablos et al. 2016a, 2016b). Thus, within the present work, MODIS Aqua daytime LST was used for estimation for 1-km daily SM in the study area.

All regression coefficients are found to be statistically significant at 5% (p < 0.05), confirming that Aqua daytime LST is a meaningful predictor for SM. In order to develop our model, the three monitoring points in the Komotini, Nimfea and Kerasia monitoring stations were used, as they cover both the plain and mountain parts and different land cover types. The three remaining points, i.e., Asomati, Kerasia, and Ifestos, were used only for the validation of the results. In that way, the developed models from the three individual monitoring points were applied for the whole study area, interpolating the regression coefficients and using IDW interpolation (RegressionCoefficients.tiff in Supplementary material). The regression equation from Equation (1) is then applied at the 1-km pixel level, computing thus daily SM estimates of the same spatial resolution for the study area.

Previous works have either used LST along with other auxiliary variables such as NDVI, soil properties, and land cover types (Mobasheri & Amani 2016; Pablos et al. 2016a; Fang et al. 2018a, 2021a) or soil temperature at the surface or at variable depths to estimate SM (Zhao & Li 2013; Ghahremanloo et al. 2019). The basic principle of these methods is that the LST is sensitive to the water content of the surface soil due to its effect on the surface heating process under the soil (Zhao & Li 2013; Mobasheri & Amani 2016). Therefore, a relationship between SM and LST can be established through the observation of soil temperature variations. (Mobasheri & Amani 2016). However, the model with LST being the only predictor variable has the ability to explain a significant part of the variability of SM (ranging from 49 to 73% according to the R² values presented in Figure 3), and there is still some variability that it is not captured by LST alone, therefore a strategy should be adopted in order to improve the model's predictive ability. Previous works showed that introducing other variables in the model, such as vegetation density and type, soil type and topography, which are found to be meaningful predictors of SM, constitutes a reasonable approach (Vergopolan et al. 2022; Ma et al. 2023). This, however, would increase model complexity in terms of data requirements. Thus, we argue that the part of the SM variability that cannot be explained by LST alone can be explained by using in situ SM observations from the existing monitoring network. An analogous approach has been presented by Ghahremanloo et al. (2019), where they integrated weather station and satellite-derived soil temperature data. In our work, we directly incorporate in situ SM data from weather stations and a residual correction approach to conclude the potential of the proposed methodology to provide reasonable daily SM estimates at 1 km.

This seems a challenging alternative as an analogous approach has been documented to work well in the mesoscale SM patterns (Ochsner et al. 2019) and thus it would be of interest to evaluate whether it can produce satisfactory results at the smaller scales, as in the case of the present work.

The results of our work indicate that the presented approach of regression and residual correction applied to estimate daily SM performs satisfactorily, according to the calculated performance metrics. Compared to other methods of SM estimation, our methodology proved to outperform both IDW interpolation and SMAP-downscaled SM. This is evident from all the metrics examined in Tables 3 and 4. A disadvantage of the present work and all approaches that use TIR remotely sensed products which are constrained by atmospheric conditions, is that the missing data in the remotely sensed products (e.g., cloud covered areas) are reflected in the produced SM. On the contrary, interpolation techniques such as IDW produce the same number of predicted SM values as the observed ones, as seen in Table 3, whereas the residual correction approach produces approximately half of the SM values compared to the interpolated and observed ones, in the three validation stations. The downscaled 1-km daily SMAP SM produces even fewer predictions. In that way, the comparison of the performance metrics among the three approaches should be examined cautiously, as the sample number is different. Furthermore, even though the SMAP based SM demonstrates the worse performance in terms of R², compared to the other two approaches, it demonstrates comparable results for the rest of the examined metrics. Considering that downscaled SMAP soil moisture does not require in situ SM observations, which is a major advantage, especially for ungauged areas, it is a useful methodology as it demonstrated good performance in an area of complex terrain and variable land cover. Its potential for improvement in applying bias correction in areas with in situ SM networks should be further investigated in future works.

Another factor that may influence the results obtained is the spacing and spatial distribution of the in situ monitoring. In the present work, the average spacing of SM monitoring points is 8 km and covers both plain and mountainous parts, different land cover and soil types; SM is estimated at the local scale at 1-km grid, while in the work of Ochsner et al. (2019) the in situ monitoring average spacing is 30 km and SM for the whole Oklahoma is estimated at the 800 m grid. In cases where the spatial distribution of monitoring stations does not cover the topography, land cover and soil variability, biases could affect the results as the regression equations will not account for a large part of the variability in SM. Additionally, our work demonstrates that the methodology described herein adequately captures the spatial SM variability, as shown in Figures 4(a)–4(c). Figure 6 also shows that all validation points performed satisfactorily in both spatial and temporal terms. The time series presented in Figure 6 shows that the simulated SM is in good agreement with the observations throughout the time periods studied. The improved description of the spatial variability of SM is also evident when comparing the results shown in Figures 7 and 8, where Figure 7 captures the spatial features inherited from the LST. On the contrary, IDW shows a smooth transition of SM across the study area, with less potential to capture the spatial patterns of SM.

Regarding the results acquired with analogous approaches elsewhere, our methodology performs with reasonable accuracy compared to previously published research in Oklahoma at the mesoscale (Ochsner et al. 2019). Other works incorporate land cover, soil type information, various vegetation indices, topographic data, in situ SM measurements and modeled SM estimates in their models and acquire high-resolution SM estimates (Ma et al. 2023). However, our work proves to be reasonably accurate while being simple in terms of data requirements, using information from a sparse network of SM sensors at the local scale. This is explained by the fact that the MODIS LST algorithm implicitly incorporates information on topography, soil types and land cover (Wan 1999, 2014). Additionally, our approach accounts indirectly for land cover, soil types and topography since it develops different regression equations for each individual SM monitoring station. Analyzing the results from the different regression formulas, higher regression slopes (stronger dependence of SM with LST) are observed for both mountain SM stations, i.e., Kerasia and Nimfea, compared to Komotini, Asomati and Karidia stations located within the plain and hilly part of the study site. The Ifestos station in the plain, although it has the highest correlation with the LST, relates to a different and shorter period. This finding demonstrates the complex role of topography and landcover on the SM and the challenges of SM disaggregation techniques, which are also highlighted in previous works (Fang et al. 2018b; Vergopolan et al. 2022; Ma et al. 2023). Despite those challenges, the approach presented herein proved to work concisely in a complex terrain, like the one examined in the Komotini study site.

A possible extension of the present work could be the retrieval of higher resolution SM using finer scale LST, e.g., Landsat LST, combined with a data fusion technique to preserve the daily time step of MODIS LST and incorporate the improved spatial resolution of Landsat missions (Hazaymeh & Hassan 2015; Li et al. 2021). A logical extension could also be the acquisition of LST from ground monitoring points if such a network operates at the fine resolution in the area of interest. Furthermore, the potential for sub-daily SM acquisitions is worth being evaluated. Although in this work daily averaged SM was correlated to daytime LST, it is possible to extend the methodology to finer temporal scales, as long as there is availability of gridded LST at improved temporal resolution.

Stationarity and seasonality of the established SM–LST relationships should also be the focus of future work, which however requires a long time series of SM information, which was not available in the present work. We believe that the regression formulas developed herein can be applied at the regional level in the wider study area, however, site-specific relationships should be evaluated in case the methodology is applied in a different area. With the evolution of climate change, however, the validity of SM–LST relationships may change so it is recommended that the regression formulas should be re-evaluated within reasonable time periods.

As mentioned above, the limitations of the present work are mainly related to cloud cover conditions which restrict many remotely sensed products. This is especially the case in high altitudes and during winter when cloud cover is more frequent. However, during summertime when SM information is mostly needed for the planning of agricultural activities, the restrictions of cloud cover are minimized. A possible extension of our approach in order to improve the LST coverage and minimize days with no LST acquisitions is to evaluate other parameters like e.g., remotely sensed precipitation, known to be related to SM (Brocca et al. 2013). However, remotely sensed precipitation is acquired at a much lower resolution, typically ∼10 km and the estimation of SM at 1 km would require downscaling of precipitation estimates prior to their use. This would increase the complexity of the approach and would introduce uncertainty in the SM estimates. Alternatively, a methodology like the one proposed by Lu et al. (2022) for gap filling combining global reanalysis data with satellite observations could be a reasonable approach to increase coverage and minimize data gaps in the produced 1-km SM estimates.

The use of other interpolation techniques should also be evaluated. In this work due to the limited number of in situ monitoring points only simple interpolation approaches were feasible. More sophisticated interpolation techniques, e.g., kriging proved to be appropriate for the description of the spatial variability of SM at the mesoscale (Ochsner et al. 2019) and its application at the local scales in areas with a higher number of monitoring stations is worth being studied. A certain limitation of the approach presented here is its dependence on in situ SM monitoring, which is not available in all areas, and in most cases, the gridded SM has to be estimated from a sparse and sometimes intermittent network. In the application presented here at the Komotini test site, the two validation points were available for a limited period of time, while the third does not cover a full year. Nevertheless, the combination of remote sensing information and a limited number of continuously operating SM sensors proved to be adequate to capture the temporal and spatial variability of SM.

The present work demonstrates a methodology to obtain gridded daily SM at the 1-km spatial resolution using information from a sparse in situ monitoring network operating in the study area. The developed approach uses the established regression equations between MODIS Aqua daytime LST and daily averaged SM for each individual in situ monitoring point, in order to interpolate them over the study area and obtain gridded daily SM at the 1 km. Furthermore, the residuals of the estimated and in situ observed SM were computed for each monitoring point and a residual correction step was applied. The obtained results were satisfactory in all performance metrics, with MAE < 0.02 m³/m³, and also better compared to SM obtained from IDW interpolation and downscaled SMAP SM. Therefore, the proposed methodology proves to be a promising alternative to the downscaling of remotely sensed soil moisture, using direct in situ SM measurements, even in the case of sparse monitoring networks, taking into account the inherent limitations of remotely sensed products. Considering that SM sensors are relatively easy to operate and have low acquisition and operational costs, the technique presented here can contribute to the estimation of the spatial variability of SM and the provision of high-resolution gridded SM products, useful for crop production and in water resources and climate assessment works.

This research was funded within the frames of the EU project titled: WATERLINE, project id CHIST-ERA-19-CES-006.

All authors contributed to the study's conception and design. Material preparation, data collection and analysis were performed by A.G., G.F., B.F., and M.K. The first draft of the manuscript was written by A.G. and V.L. and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.