A simulation of soil water content based on remote sensing in a semi-arid Mediterranean agricultural landscape

This paper shows the application of a water balance based on remote sensing that integrated a Landsat 5 series from 2009 in an area of 1,300 km 2 in the Duero Basin (Spain). The objective was to simulate the daily soil water content (SWC), actual evapotranspiration, deep percolation and irrigation rates. The accuracy of the application is tested in a semi-arid Mediterranean agricultural landscape with crops over natural conditions. The results of the simulated SWC were compared against 19 in situ stations of the Soil Moisture Measurement Stations Network (REMEDHUS), in order to check the feasibility and accuracy of the application. The theoretical basis of the application was the FAO56 calculation assisted by remotely sensed imagery. The basal crop coefficient (K cb ), as well as other parameters of the calculation came from the remote reflectance of the images. This approach was implemented in the computerized tool HIDROMORE+, which integrates various spatial databases. The comparison of simulated and observed values (at different depths and different land uses) showed a good global agreement for the area ( R 2 = 0.92, RMSE = 0.031 m 3 m –3 , and bias = –0.027 m 3 m –3 ). The land uses better described were rainfed cereals ( R 2 = 0.86, RMSE = 0.030 m 3 m –3 , and bias = –0.025 m 3 m –3 ) and vineyards ( R 2 = 0.86, RMSE = 0.016 m 3 m –3 , and bias = –0.013 m 3 m –3 ). In general, an underestimation of the soil water content is noticed, more pronounced into the root zone than at surface layer. The final aim was to convert the application into a hydrological tool available for agricultural water management.

and spatially distributed (Zhang & Wegehenkel, 2006;Er-Raki et al., 2010;Sánchez et al., 2010). The use of a linear relationship of NDVI-K cb (Normalized Difference Vegetation Index-basal crop coefficient) instead of the tabulated values of K cb improves the accuracy of the results, even over woody crops, as shown in Campos et al. (2010). Another promising aspect for using the NDVIbased K cb is its ability to provide a direct approach to determining actual evapotranspiration (AET) conditions when crop growth and water use deviate from optimum conditions (Hunsaker et al., 2005).
The second basic approach in the application is using a map of land use/land cover extracted from an image classification method. The map can produce spatial distributions of some calculation parameters assigned to each of its classes.
In areas where irrigation rates must be controlled precisely, where it is essential to prevent water loss, one approach to check if the supply is adequate to meet the demand is to account for the respective components in the water balance (Singh et al., 2010). In this sense, irrigation advisory systems show increasing success in different regions of Spain and the world based on climate information, crop status, and agricultural practices (Martín de Santa Olalla et al., 2003).
This work aims to demonstrate the reliability of the combination of remote sensing with a hydrologic balance model to obtain hydrological variables in a distributed manner. For this proposal, this work compares the results of SWC simulated against in situ measured values. A second objective is to demonstrate the accuracy of this straightforward method so that it can be implemented as a query tool for farmers and managers of water, who in the future could use it to obtain near real-time data of crops and water requirements. Here, the water balance is implemented in a spatially distributed tool, HIDRO-MORE+, which was previously tested in the same area in 2001 and 2002 (Sánchez et al., 2010).

Introduction
Soil water content (SWC) is a highly dynamic ecological variable, but it is also one of the most sensitive factors in agricultural yield. In semi-arid Mediterranean areas, where water scarcity is a common situation, SWC is the limiting factor for agriculture, and crop evapotranspiration is the main factor of water consumption. Accurate information on SWC is crucial for practical agricultural water management at various scales (Ju et al., 2010). Point measurements of SWC are frequent in agricultural areas for controlling irrigation scheduling and water reserves (Bonet et al., 2010), either with measurements at different depths or with a water balance using the crop coefficient-reference evapotranspiration (K c -ET 0 ) approach (Allen et al., 1998). However, the spatial heterogeneity of soil properties makes it difficult to scale up measurements from points to large scales, and spatial monitoring designs are costly and timeconsuming (Western & Grayson, 2000). This is a difficult problem for agricultural water management, and several remote sensing spatial missions (such as AMSR-E, ERS, SMOS, and SMAP) have attempted to fill this gap, even though improved process understanding and algorithms are needed to enable the use of these data in the future (Wagner et al., 2007a).
The application of hydrological models to simulate SWC with adequate accuracy in heterogeneous and complex areas remains challenging, and the utmost efforts should be made to exactly determine fine-scale meteorological data and to build a spatial skeleton to improve soil moisture modeling (Rößler & Löffler, 2010). The possibility of producing maps of hydrological variables by combining a mass balance approach and remotely sensed data to calculate a water budget is a useful tool for watershed-level hydrologic models (Earls et al., 2006) because the water balance provides high temporal resolution, whereas the remote sensing approach provides high spatial coverage. It is thus possible to monitor crop status and water use through remotely sensed crop coefficients in a dynamic manner for every plot. The crop coefficients derived from the vegetation index (Bausch & Neale, 1987;Neale et al., 1989;Calera et al., 2004;Gonzalez-Dugo & Mateos, 2008;Er-Raki et al., 2010) account for any situation out of pristine growing conditions of climate, management or vegetative vigor.
The FAO56 water balance (Allen et al., 1998) has been used in recent works both at point scale (Allen, 2000;Hunsaker et al., 2005;Er-Raki et al., 2007; A simulation of soil water content based on remote sensing automatic weather stations among other equipment for hydrological and climatic monitoring (Martínez-Fernández & Ceballos, 2003. These stations are located in an area of approximately 1,300 km 2 (41.1° -41.5° N; -5.1° -5.7° W) in a central semi-arid part of the Duero Basin (Fig. 1). This area is nearly flat, ranging from 700 to 900 m a.s.l. The climate is continental semiarid Mediterranean. The land uses are mainly agricultural: rainfed cereals in winter and spring (78%), irrigated crops in summer (5%), and vineyards (3%). There are also some patchy areas of mixed forests and pastures (13%). The most abundant soils are Luvisols and Cambisols with a predominantly sandy texture (mean sand content, 71%) (Martínez-Fernández & Ceballos, 2003).
In the year of study (2009), the rainfall was 332.1 mm, and the ET 0 (FAO56-Penmann-Monteith method) was 1,147.6 mm. The monthly data (Fig. 2) show a water deficit for all months except December and January. ET 0 is especially high in the summer period.

Input data and pre-treatment
The simulation was performed along 2009 at daily interval. Three datasets were needed for the calculation: meteorological information (i.e., daily precipitation and ET 0 ), imagery (i.e., NDVI time series and a map of land use-land cover), and spatially distributed soil information.
Daily precipitation and ET 0 , were computed from the automatic weather stations (Fig. 1). HIDRO-MORE+ performs a spatial interpolation based on the Inverse Distance Weighting method (Shepard, 1968).
The imagery came from Landsat 5 TM (scene 202-031). In 2009, three images were selected after discarding cloudy images: March 23, June 11, and August 30. These dates covered the most significant periods of the crop growing cycles under study, i.e., spring for rainfed cereals and summer for irrigated crops and vineyards. The images were orthorectified and registered by a rigorous model (Toutin, 2004) with 20 ground control points, plus the elevation digital model and the orbital information. This process guarantees the spatial coincidence of the three scenes. Radiometric calibration and atmospheric correction were also performed by means of the standard values of the atmosphere and the dark object subtraction method (Chavez, 1989). The NDVI (Rouse et al., 1974) was calculated as the normalized difference between infrared and red physical units, i.e., the  reflectance at-surface resulting from the radiometric and atmospheric treatment. The composite image of the three NDVI become the input for the classification process and the basis for computing K cb as well as the fraction of vegetation cover and other parameters of the calculation. For the classification, a hybrid method was used that mixed the Geographic Information System (GIS) for Agricultural Plots of Castilla y León Region (SIGPAC) with the NDVI composite. This process comprised two steps. First, the SIGPAC plots of water, urban, forestpasture and vineyards were clipped from the image. Second, to separate the mixed category 'arable land', a segmentation procedure of the NDVI levels was performed (Vincent & Pierre, 2003), which allowed the segmentation of rainfed, irrigated and fallow areas, taking advantage of the opposite growing cycle of all of them. The union of the results of both steps became the land use/land cover map (Fig. 3). The accuracy assessment with ground truth areas showed a coincidence of 100% for water, 85.7% for forest-pasture, 76.5% for vineyard, 84.6% for rainfed cereals, 100% for irrigated crops and 90% for fallow. The average accuracy was 86.8%. This map established some parameters of the water balance for each class (or restrict some others), such as root depth, plant height, and thresholds of K cb .
Finally, for the soil database required for the calculation, soil samples at surface level (0-5 cm) were collected in a grid basis (cells of 3 km × 3 km). The SWC at field capacity (SWC FC ) and at wilting point (SWC WP ) were determined for the center of 146 cells. The water retention curve was obtained in laboratory measurements on these undisturbed soil samples (100 cm 3 monoliths) taken at the surface layer in the center of each cell. The method combines tension boxes (sand box method) for potentials close to zero and membrane pressure method for low tension conditions. Using this combined method nine experimental points (soil moisture vs. tension pairs) of the soil retention curve were obtained at 0, 0.2, 1, 3, 10, 30, 100, 300, and 1,500 kP. After the adjusting to an empirical model (van Genuchten, 1980) one can obtain the entire curve and then obtain parameters as field capacity (-33 kPa for SWC FC ), wilting point (-1,500 kPa for SWC WP ), saturation or so on.

Calculation basis
The K c -ET 0 approach uses meteorological data and crop coefficients to calculate crop evapotranspiration. The dual form of K c was chosen (Wright, 1982), as is developed in FAO56: The term ET 0 K s K cb represents the transpiration component from the plants, and the term ET 0 K e , represents the evaporation component from the soil. K cb is the transpiration coefficient at a potential rate, i.e., when water is not limiting transpiration, and it is usually obtained from tabulated values; here, it was calculated as a linear form of the NDVI (Bausch & Neale, 1987): K s describes the effect of water stress, and it is calculated according the water content in the root layer, estimated from the daily water balance. The soil evaporation coefficient, K e , is calculated from the daily water balance on the upper soil surface layer (10 cm). Thus, the calculation accounts for the effective root zone, but the upper topsoil is considered during the non-growing periods or over bare soil conditions (i.e., fallow). This daily balance, the coefficients calculation, and the limits of available soil water content are described by Allen et al. (1998) and Allen (2000). SWC can be calculated as the residual of the daily water balance. It is expressed in volumetric units and is analyzed in terms of plant available water. FAO56 establishes the SWC FC as the upper limit, and the remaining SWC is calculated as the difference between the water at field capacity and the water deficit or depletion (the residual of the balance). HIDROMORE+ establishes the SWC WP as the minimum SWC.

Comparison methods and alternatives
Hydrological models are mainly calibrated and validated using runoff (Rößler & Löffler, 2010) or evapotranspiration (Allen, 2000;Eitzinger et al., 2002;Er-Raki et al., 2007). However, because runoff processes do not provide much insight into internal processes (Grayson et al., 1992a,b) soil moisture simulations also need to be validated. Regarding this validation purpose, the REMEDHUS network performs a continuous measurement of 0-5 cm soil moisture in 19 stations with one Hydra Probe (Ste-vens® Water Monitoring System, Inc.), and at different depths with an EnviroSMART TM soil water profile (Sentek Technologies) installed in 12 stations. The EnviroSMART probes can include several sensors that measure soil water at multiple depths, customized for individual applications by adding sensors at user-specified measurement. In REMEDHUS, three probes at 25, 50 and 100 cm depths were installed. Hydra and EnviroSMART are capacitance probes which measure soil dielectric constant to calculate soil moisture content. At each station sensors are connected to a CR200 (Campbell Scientific Ltd.) datalogger where soil moisture data are stored at hourly intervals. As stated before, there were 19 Hydra stations in that period; among them, 12 also had EnviroSMART profiles.
The hydraulic properties were obtained in the laboratory in a similar procedure as the soil grid, i.e., after taking undisturbed soil samples at the different depths in each station. SWC FC , SWC WP and SWC at saturation (SWC sat ) were determined. For the validation process, taking into account that water balance considers the SWC at root depth, profile-average values were considered as well as the root depth for each land use/land cover (Table 1).
The alternatives considered for the SWC validation were the following: a) using the surface SWC observed in the Hydraprobe stations (SWC hyd ); b) using the profile-averaged SWC resulting from the Enviro-SMART profiles (SWC env ); c) using the profile-averaged SWC of all sensors (one at surface and three along the profile), SWC all .
Six stations under rainfed cereals, five under vineyards, five under fallow and three under forest-pasture were used. The validation strategy consisted in comparing simulated and observed values using the coefficient of determination (R 2 ), bias, and the root mean square error (RMSE), namely:

bias SWC SWC
where the bias represents the deviation or difference between the simulated SWC and the observed series of SWC (SWC hyd , SWC env , and SWC all ) for each day i; and the RMSE represents the cuadratic mean of these differences along the time series. These statistics were applied for each station, as well as for the land use-averaged SWC and total-averaged SWC of the whole area.

Results
The simulated SWC was compared against the observed following the three alternatives of validation (Fig. 4, only one station representative of each category is shown). The vineyard soils took small SWC values along the year (Fig. 4a). These soils have more content of sand and thus less retention capacity, and there are no differences among the observed values along the profile. The simulated SWC curves and the observed curves matched perfectly. Therefore, for the other categories, there are greater differences between SWC hyd , SWC env , and SWC all . Among these differences, as expected, the behavior of SWC hyd fluctuates more because it is a surface layer content; this fact is more evident in the rainfed cereal (Fig. 4b). The failure of the simulation of fallow lands is noticeable (Fig. 4c). Because the calculation does not consider plant activity, but only the evaporation layer, the depletion is the maximum and SWC equals the minimum value, SWC WP .
For the simulated SWC, a slight underestimation compared with the observed values can be seen in all categories, especially if comparing SWC hyd with the forest-pasture cover (Fig. 4d). The observed SWC is greater in forest-pastures due to their location in the valley bottoms, where the water table is shallow in winter, floods sometimes occur, and there is a higher content of clay in the soils. It should be pointed out that the calculation establishes the upper limit of the simulated SWC as SWC FC , but the curve indicates that the actual SWC on the soils exceeds this limit in at least two situations: a) among the deep profile values (SW-C all and SWC env ) and b) among the surface values (SWC hyd ) during the rainy periods.
The results of the comparison (Fig. 5) were considered satisfactory in terms of coefficient of determination (R 2 > 0.6) for all stations except for the fallow category and one isolated station (I6). There are no meaningful differences between the three alternatives of validation in terms of determination coefficient (Fig. 5a). Therefore, we can conclude that the simulation is useful for correctly characterizing the mean profile SWC as well as the surface level SWC. But the higher RMSE at Fig. 5b indicated a worse fitting with the surface values in most of the stations, and the fluctuating sign of the bias at this level (Fig. 5c) indicates that the simulation may barely predict the higher variability at the surface, whereas the bias for the profile-averaged content is always a dry bias.
The estimation improved if the simulated SWC was averaged for each land use/land cover and for the whole area ( Table 2). The area-averaged simulation agrees bet- ter with both the mean observations at surface level (R 2 = 0.94, RMSE = 0.041 m 3 m -3 , and bias = -0.002 m 3 m -3 ) and the profile-average (R 2 = 0.86, RMSE = 0.037 m 3 m -3 , and bias = -0.034 m 3 m -3 for SWC env ; and R 2 = 0.92, RMSE = 0.031 m 3 m -3 , and bias = -0.027 m 3 m -3 for SWC all ). As stated before, the underestimation is still clearer for the profile-averaged than for the surface values. Regarding the categories best characterized, the simulation is slightly better for vineyards and for rainfed cereals. The comparison cannot be done over forest-pasture because only one EnviroSMART profile is located in this land use (H9), which represents an isolated value. However, the error for this category at surface level is the largest (RMSE = 0.117 m 3 m -3 ).

Discussion
Many researchers have studied the relationship between root growth and soil water (Quanqi et al., 2010) because uncertainties come mainly from the high variability of soil properties and root depth (Ju et al., 2010;Sánchez et al., 2010). These parameters of soil properties and root activity have shown to be the most important parameters affecting the SWC on a small scale. The dependence of root activity in this model is verified here through the failure of the simulation over fallow or bare soil plots. The present work used SWC as control variable of the water balance in a profilebased observation. Since the in situ observations were made at different depths, the single value of simulated SWC can be checked against the vertical profile; in particular we used the upper topsoil values and the profile-averaged. The simulation represents fairly well the SWC dynamics compared with both the surface curve and the soil profile, but the bias is reduced when comparing the simulated SWC with the profile-averaged SWC (SWC all ). The underestimation observed in Sánchez et al. (2010) was also detected in this case because the profile-averaged SWC is clearly underestimated in eight of the twelve stations. As suggested there, there is a need to redefine the limits of the plant available water used in the calculation, which in turn depends on the field capacity and the wilting point. However, the RMSE obtained here, ranging between 0.011 and 0.149 m 3 m -3 , is in line with similar works (Diekküger et al., 1995;Vanclooster & Boesten, 2000;Jalota & Arora, 2002;Zhang & Wegehenkel, 2006), and with remotely sensed estimations of soil moisture, such as SMOS (Soil Moisture and Ocean Salinity from the European Space Agency) and SMAP (Soil Moisture Active-Passive from NASA), targeted in an RMSE of 0.04 m 3 m -3 . It is necessary to note the higher RMSE values in the case of the three forest-pasture stations at surface layer. The RMSE was 0.149, 0.085 and 0.164 for H9, K9 and M13, respectively. These errors can be explained for the specific location of these stations in the bottoms of the valleys, near the stream water, where the water table is shallow in winter, flooding occasionally occurs, and there is a higher content of clay in soils. In these cases, the actual SWC could reach higher values than the SWC FC in some periods, and these conditions cannot be explained by the standard calculation of FAO56. Another explanation is that the FAO56 method is generally limited to applications with crops, although it can be applied to natural vegetation with greater uncertainty in the crop coefficient (Allen & Bastiaanssen, 2005) or woody crops such as vineyards (Campos et al., 2010). However, experience suggests that the most influential parameter in the water balance is the water retention capacity, expressed in terms of SWC FC , and the root depth. In the application of a distributed, physically based model, Rößler & Löffler (2010) found that soil moisture was especially sensitive to changes of soil parameters, and here, the most critical parameter seemed to be the SWC FC , the upper limit of plant water availability.
In previous experiments of soil moisture estimations in REMEDHUS Wagner et al., 2007b;Sánchez et al., 2010), the soil moisture was observed every two weeks, whereas a daily averaged interval was used here. This new interval matches the time interval of the calculation and improves the reliability of the validation process because the potential variability of both curves is expected to be similar. In the present study, the simulated soil moisture shows a good agreement (total-averaged SWC with R 2 = 0.94, 0.86 and 0.92 for SWC hyd , SWC env and SWC all , respectively), whereas the values are mostly underestimated (bias negative), especially in the . This underestimation is in general clearer when comparing the simulated SWC with the whole profile average due to its higher water content. For the same reason, the underestimation was more pronounced during the wet season, i.e., winter; reaching bias above -0.40 m 3 m -3 for forested stations such as H9 in this period (Fig. 4). High values of water content may have been hindered by the limit of the SWC FC established in the simulation. It can be concluded that the simulation better described the dynamics of surface SWC, in terms of R 2 , than at the deeper layers, but the errors are higher (Table 2) and more variable (Fig. 5). The best simulation took place for vineyards and rainfed cereals. One limitation is that the relationship between remote VI and crop coefficients were based on the basal crop coefficient K cb , which is defined when the soil surface is dry but transpiration is occurring at potential rate (Allen et al., 1998). These conditions are particularly present in herbaceous crops (e.g., cereals), due to this high fraction of coverage (minimum soil evaporation), and during the growing cycle (i.e., March to June), free of water stress in the root zone. This fact could explain the good results of the approach in this area, due to the predominance of this coverage, and the poor agreement over the fallow areas. These con ditions are not present over vineyards, on the contrary, there is a wide area of bare soil among plants, with scarce water supply during their development stage. Even though, the simulation seems to be quite accurate for these stations in average (R 2 = 0.92, RMSE = 0.031 m 3 m -3 , bias = -0.027 m 3 m -3 ). It could be possible that the low water contents in these sandy soils would mask the errors of the estimation. Working under non-stressed and reduced evaporation conditions, Campos et al. (2010) obtained a good agreement of simulated AET over vineyards using a linear relationship NDVI-K cb . These findings support the possibility of an accurate application of the remotely sensed K c -ET 0 methodology, and it opens the suite of crops for which its application is suitable.
The hypothesis of the suitability of HIDROMORE+ and the NDVI-K cb approach for the distributed calculation of hydrological variables, particularly SWC, has been tested using the database of the REMEDHUS network (Duero basin, Spain). The results have been very satisfactory, especially in agricultural areas (vineyard and rainfed cereals), in terms of averaged SWC for the whole area. In the fallow areas, where the root activity is absent or limited, the calculation failed. Indeed, the critical parameters of the SWC calculation have proven to be the root depth and the upper limit of the SWC considered by the model. This limit, established as the water content at field capacity, prevents higher values of SWC in the simulation. Hence, the forest-pasture cover showed an underestimation of the surface SWC, and a general underestimation around -0.030 m 3 m -3 is detected for all the stations in the root zone, more clear during the wet period.
The procedure described in this paper for estimating soil water content and other components of the water balance opens the possibility of combining the measured point-data with imagery products from other optical, thermal or radar bands and sensors to increase the accuracy of the spatially distributed water balance. Obtaining an accurate distributed water balance is the main hydrological challenge for many practical and theoretical purposes related with water management and irrigation schemes along a wide range of scales. Table 2. Determination coefficient, R 2 , root mean square error (RMSE, m 3 m -3 ), and bias (m 3 m -3 ) for each category (averaged and the total) averaged in the three alternatives of validation. There is no average for forest-pasture because there is only one station with both sensors their support (ESP2007-65667-C04-04, AYA2010-22062-C05-02 and EBHECGL2008-04047 projects), which made this work possible.