Analyzing Spatio-Temporal Factors to Estimate the Response Time between SMOS and In-Situ Soil Moisture at Different Depths

A comprehensive understanding of temporal variability of subsurface soil moisture (SM) is paramount in hydrological and agricultural applications such as rainfed farming and irrigation. Since the SMOS (Soil Moisture and Ocean Salinity) mission was launched in 2009, globally available satellite SM retrievals have been used to investigate SM dynamics, based on the fact that useful information about subsurface SM is contained in their time series. SM along the depth profile is influenced by atmospheric forcing and local SM properties. Until now, subsurface SM was estimated by weighting preceding information of remotely sensed surface SM time series according to an optimized depth-specific characteristic time length. However, especially in regions with extreme SM conditions, the response time is supposed to be seasonally variable and depends on related processes occurring at different timescales. Aim of this study was to quantify the response time by means of the time lag between the trend series of satellite and in-situ SM observations using a Dynamic Time Warping (DTW) technique. DTW was applied to the SMOS satellite SM L4 product at 1 km resolution developed by the Barcelona Expert Center (BEC), and in-situ near-surface and root-zone SM of four representative stations at multiple depths, located in the Soil Moisture Measurements Station Network of the University of Salamanca (REMEDHUS) in Western Spain. DTW was customized to control the rate of accumulation and reduction of time lag during wetting and drying conditions and to consider the onset dates of pronounced precipitation events to increase sensitivity to prominent features of the input series. The temporal variability of climate factors in combination with crop growing seasons were used to indicate prevailing SM-related processes. Hereby, a comparison of long-term precipitation recordings and estimations of potential evapotranspiration (PET) allowed us to estimate SM seasons. The spatial heterogeneity of land use was analyzed by means of high-resolution images of Normalized Difference Vegetation Index (NDVI) from Sentinel-2 to provide information about the level of spatial representativeness of SMOS observations to each in-situ station. Results of the spatio-temporal analysis of the study were then evaluated to understand seasonally and spatially changing patterns in time lag. The time lag evolution describes a variable characteristic time length by considering the relevant processes which link SMOS and in-situ SM observation, which is an important step to accurately infer subsurface SM from satellite time series. At a further stage, Remote Sens. 2020, 12, 2614; doi:10.3390/rs12162614 www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 2614 2 of 28 the approach needs to be applied to different SM networks to understand the seasonal, climateand site-specific characteristic behaviour of time lag and to decide, whether general conclusions can be drawn.


Introduction
Soil moisture (SM) is recognized as an essential climate variable (ECV) and it plays an important role in agricultural applications including drought monitoring and rainfall estimations, weather forecasting and climate models across different spatial and temporal scales. Therefore, comprehensive knowledge about SM content has substantial socio-economic relevance to food security. Meteorological conditions set the boundary conditions at the atmosphere-soil interface for near-surface SM [1]. Most of cropping systems develop their largest amount of root biomass within the first 50-100 cm depth of the soil, which is commonly defined as the root-zone [2]. Surface water balance comprises processes such as infiltration, runoff, evaporation and transpiration, and subsurface water movements including recharge and drainage. These processes occur at different temporal scales and knowledge about their variability is critical to effectively determine crop growing seasons and to schedule irrigation if necessary [3]. For increasing climate variability, long-lasting events such as heat waves and precipitation hold off become more frequent [4]. An understanding of both near-surface and root-zone SM is of key importance to study water scarcity and drought phenomena [5,6].
SM acquisition methods can be generally divided into air-and spaceborne remote sensing techniques, non-invasive geophysical measurements such as Ground Penetrating Radar (GPR) and Electromagnetic Induction (EMI), and invasive in-situ observations [7][8][9][10]. Spaceborne remote sensing techniques for SM retrieval are relatively inexpensive and provide high spatial coverage which cannot be achieved using labor-intensive geophysical measurements [11]. Besides that, in-situ observations enable the continuous monitoring of SM from multi-depth measurements at point-scale. Therefore, the use of satellite imagery in combination with more accurate point-scale observations seems to be suitable for analyzing SM variability. The Soil Moisture and Ocean Salinity (SMOS) mission, launched in 2009 by the European Space Agency (ESA), is the first space mission specifically dedicated to estimate SM. It uses an L-band interferometric radiometer to provide global SM with a revisit time at the equator of three days and a spatial resolution of ∼40 km [12]. The Soil Moisture Active Passive (SMAP) is the second spaceborne mission devoted to monitor SM [13]. It was launched in 2015 by the National Aeronautics and Space Administration (NASA) and carries on board an L-band real aperture radiometer. It was also equipped with a radar, but it failed after three months in orbit. SMAP and SMOS provide global SM at ∼40 km every three days. Among the several in-situ SM networks available worldwide and included in the International Soil Moisture Network (ISMN) [14], the REMEDHUS network of the University of Salamanca in Western Spain provides SM observations at multiple depths with excellent temporal resolution [15,16]. This network has been widely used for calibration and validation of both passive and active remote sensing SM data [17][18][19][20][21].
Remote sensing measurements are sensitive to SM content within a certain layer and over a particular area. The depth range of this layer within the soil profile is defined by the penetration depth and the extension of the area is determined by the sensor footprint size. Thus, the validation of satellite measurements requires several assumptions. Firstly, frequencies at L-band (1-2 GHz) are sensitive to SM at approximately the top 5 cm of the soil [22], commonly known as surface SM. However, the penetration depth does not only depend on the frequency of the electromagnetic signal, decreasing as frequency increases, but also on the attenuation of the signal due to changes in the soil temperature and the SM. It is larger in dry than in wet soils [23]. Secondly, satellite SM is usually compared against in-situ observations acquired by sensors installed at 5 cm depth, using collocated and concurrent data. Notwithstanding, both SMOS and SMAP SM respond more quickly to wetting up and dry down events, similar to the precipitation pattern, revealing a response time between satellite and in-situ measurements at the top 5 cm [24][25][26][27]. Finally, satellite SM measurements are assumed to have low intra-pixel variability. To deal with the issue of comparing data at different spatial scales, different strategies have been used to upscale in-situ SM [28][29][30]. A further assumption in validation is that the location of the point-scale SM observation is presumed to be most representative for the surrounding area, which is considered homogeneous.
The unsaturated (or vadose) zone is the soil layer between the land surface and the groundwater table, including the capillarity fringe, and its thickness varies from zero meters in the lakes to hundreds of meters in arid regions. The hydraulic conductivity, which describes the water movement through the soil porous media, depends on soil properties including porosity, permeability and saturation. Except in given cases of water uptake, SM tends to move mainly downwards due to gravity. The rate of SM exchange can be defined by two separable fluxes corresponding to advective and diffusive terms [31]. The diffusive term describes the exchange of SM due to the shape of the soil water capillarity profile and becomes important under long-lasting dry conditions, where hydraulic conductivity is low. In case a sharp wetting front enters the soil after heavy rains, the diffusive term becomes negligible and the advective SM exchange dominates during water infiltration. Estimation of SM along the soil depth profile requires knowledge about natural SM-related processes, such as precipitation, evapotranspiration, infiltration and surface runoff, soil layering and man-made processes, such as irrigation. For instance, the distribution of SM and soil properties determine how much precipitation infiltrates into the soil or results in runoff. These two related processes affect the near-surface soil layer at relatively large spatial scale and at different timescales [32,33]. Additionally, small-scale variations in topography and soil properties (local heterogeneity) are important factors regarding the local distribution of SM, which makes surface SM highly variable. Consequently, surface SM observations are mainly characterized by short-term fluctuations, while root-zone SM is less dynamic, more affected by long term atmospheric conditions and dominated by seasonal trends.
Continuously changing atmospheric conditions and corresponding processes are supposed to be reflected in the time series of satellite observations. Nowadays, satellite SM retrievals are widely used to monitor surface SM [34][35][36], but they may also contain useful information about the temporal SM variability of the underlaying soil profile to infer root-zone SM within a preceding time span. Based on this, a simple method to estimate subsurface SM trend was developed by using an exponential filter, which was applied to remotely-sensed surface SM time series without the availability of local soil properties and complete drainage information [37]. Hereby, subsurface SM could be quantified by the Soil Water Index (SWI), normalized between 0 and 1, representing the relative subsurface SM dynamics at different depth levels. The behavior of the SWI is mainly described by a characteristic time length (also known as T parameter). It indicates the degree of smoothing applied to surface SM and depends on the selected depth and a pseudo-diffusivity constant controlling the flow between the soil layers.
This time length is considered as a surrogate parameter which encompasses all relevant processes and soil properties that influence the SM variability at different temporal scales. Root-zone SM was estimated from different active and passive satellite SM retrievals (including European Remote-Sensing Satellite (ERS)-Scatterometer, Advanced Scatterometer (ASCAT), Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E) and SMOS) and the obtained SWI was evaluated at several in-situ SM networks over Europe [18,[37][38][39][40][41]. Hereby, the characteristic time length was optimized as being constant over the course of a year or even some years for specific location and depth level. However, the response time in which subsurface SM is affected by atmospheric forcing varies, being shorter around precipitation events and longer during dry events [42]. Consequently, the characteristic time length cannot be assumed as constant, but seasonally variable. This variability could be particularly notable if the water balance becomes perturbed in case of high rate of water introduction or removal due to pronounced precipitation (or irrigation), or strong evaporation and transpiration.
There are different approaches to quantify the time lag between non-linearly related time series, which in principle can be divided into representation methods and similarity methods [43]. In the first case, shorter sub-sequences can be approximated from longer time series with the aim to identify representative patterns which relate the fundamental features of different observations [44,45]. Wavelet transforms can be analyzed to decompose time series into similarly shaped sub-sequences. As an example in hydrology, singular value decomposition (SVD) was applied to identify complex spatial response modes of temporal sequences, and phase information of the cross-wavelet was evaluated to study the response of groundwater levels to precipitation [46]. Their time series were related to extract anomalies and extreme events such as heavy rainfall, and local time lags were effectively approximated. In contrast, similarity methods aim for sequence alignment based on a distance measure which describe the point-to-point alignment between two time series, rather than relating single pronounced events or representative sequences. Dynamic Time Warping (DTW) is a method capable of evaluating the similarity between two time series. It analyses local variations, distortions and shifts between the time series and quantifies the time lag by finding the optimal match between them. DTW and derived versions provide the most effective similarity measure to align time series data in a wide variety of applications [43,47]. The temporal evolution of the time lag may provide information about mutual dependencies or causal relationships between the time series. Originally, DTW was applied in speech recognition to identify a particular word within a longer, more distorted signal [48,49]. In the remote sensing community, DTW was mostly used for land cover classification purposes [50][51][52]. In addition, it was used to analyze similarity of different satellite SM estimations to in-situ near-surface SM observations [53].
In this study, the DTW technique was customized and applied to estimate the response time between remotely sensed and in situ SM at different depths. The objective was to investigate and quantify the relative temporal behavior of these time series, which can be non-linearly related, to better understand the link between satellite and subsurface SM. The SMOS Level 4 (L4) SM data product at 1 km resolution provided by the Barcelona Expert Center (BEC) [54] and in-situ SM measured at several depths over REMEDHUS [19] were used. Since the temporal resolution of SMOS is limited by its revisit time, time series were averaged over three days. They were normalized to account for differences in dynamic ranges between area-averaged satellite and point-scale in-situ observations. The observations are related based on the most prominent common features contained in both time series.
Satellite time series were compared to SM observations at 5, 25 and 50 cm depth at four representative in-situ stations to estimate the time lag between surface and near-surface as well as root-zone SM, respectively. Most of the REMEDHUS area consists of cultivated fields and the interaction between SM and vegetation is non-negligible regarding the influence of crops during their growing season [55,56]. Therefore, spatial heterogeneity of land use was analyzed to provide information about the level of representativeness of SMOS observations to each in-situ station. Land use was investigated by means of the Normalized Difference Vegetation Index (NDVI) using composite time series data from the Moderate Resolution Imaging Spectroradiometer (MODIS) and high-resolution images from Sentinel-2 L1C [57]. To support the interpretation of the time lag evolution, SM seasons were estimated from precipitation and evapotranspiration recordings to categorize sub-periods in which similar SM-related processes occur. The evolution of the quantified response time was analyzed considering the results of the spatio-temporal analysis of the study area to recognize seasonal and spatial patterns, and to give insight into the variability of the characteristic time length.
The SMOS mission official L3 and the high-resolution SMOS L4 SM product are described in Section 2.1.1. The REMEDHUS SM network including multi-depth observations under consideration of the predominant climate factors are described in Section 2.1.2. DTW is introduced in Section 2.2; its fundamentals are explained in Section 2.2.1. Customization of the method including warping path step-size condition, maximum allowed time lag and the determination of onset dates of pronounced precipitation events are discussed in Section 2.2.2. Results of the estimation of SM seasons as part of the analysis of temporal variability of climate factors and the evaluation of the spatial heterogeneity of land use are given in Sections 3.1 and 3.2, respectively. Intermediate steps of the application of DTW and customization to SMOS and in-situ time series are assessed and presented in Section 3.3. The resulting evolution of time lag is presented and interpreted in Section 3.4. The potential of the DTW technique regarding the interpretation of SM processes and the inference of subsurface are discussed in Section 4. Major findings are summarized in Section 5.

SMOS Soil Moisture
There are several ongoing spaceborne missions for SM retrieval using active (radar) and passive (radiometry) sensors at microwave frequency including ASCAT (5.255 GHz at C-band), SMAP and SMOS (∼1.4 GHz at L-band), and the Advanced Microwave Scanning Radiometer 2 (AMSR-2) (frequencies ranging between 6.9 to 89 GHz), providing SM products with a spatial resolution of several tens of kilometers at global coverage [58]. ESA's SMOS mission was launched in November 2009 and is the first spaceborne mission specifically designed to provide global SM and sea surface salinity maps [59,60]. SMOS is equipped with the Microwave Imaging Radiometer with Aperture Synthesis (MIRAS), an interferometric radiometer that acquires multi-angular full polarization brightness temperatures in ascending and descending orbits with equator crossings at 06:00 and 18:00 local times [61]. It has a nominal accuracy of 0.04 m 3 m −3 and a spatial resolution of ∼40 km. A 2-3 days revisit time is sufficient to characterize the root-zone SM in case auxiliary meteorological information is included [59]. The SMOS retrieval algorithm is based on an iterative approach that minimizes a cost function by finding an optimum set of SM and vegetation characteristics using multiple incidence angle brightness temperatures [62]. Hereby, the extraction of SM from brightness temperature consists of the following steps. The observed brightness temperature is normalized to emissivity using a radiative transfer model to account for effective surface temperature and the effects of soil surface roughness. A tau-omega model is used to correct for the attenuation of soil emission due to a vegetation layer. The emissivity measurement is related to the soil dielectric properties, which are changing significantly under variation of water content regarding the medium and are estimated at any surface [62]. In a final step, the corresponding properties are related to SM.
For the last 10 years, several products with enhanced spatial resolution have been developed [63]. In this study, time series of the high-resolution SMOS L4 SM product at 1 km provided by Barcelona Expert Center (BEC) have been used [54]. It is obtained by a linear downscaling algorithm-similar to Reference [64]-merging data from microwave and visible infrared sensors with different spatial resolutions, and modeled data. The relationship between SMOS L3 SM product and auxiliary data is derived from the so-called universal triangle method in the vegetation-land-surface temperature space, using an adaptive moving window to ensure homogeneous climate conditions. The auxiliary data comprise SMOS L1C brightness temperatures, NDVI from MODIS, and land surface temperature of the European Center for Medium-Range Weather Forecasts (ECMWF). Previous research revealed a slight difference between SMOS brightness temperatures measured at ascending and descending orbits over the same homogeneous scene [65], resulting in non-negligible differences in SMOS SM retrievals. Therefore, both orbits are only compatible for specific applications. In the following, high-resolution SM maps of only the ascending orbits were used.

REMEDHUS Soil Moisture and Climate Data
The REMEDHUS network is located in the Duero basin comprising a relatively flat area of 35 km × 35 km between 41.1-41.5 • N and 5.1-5.7 • W ( Figure 1) with small slopes and an altitude ranging between 700-900 m above sea level. Land use has remained broadly stable along the last decade and it consists predominantly of rainfed cereals (80%) (winter wheat and barley), forest pasture (12%), irrigated crops (5%) and vineyards (3%) with occasional legumes and sugar beet. The growing season of winter wheat is usually between October and July, but can vary by several weeks depending on the meteorological conditions. The main soil texture consists of sandy loam and sandy clay loam with the topsoil becoming sandier at the flat spurs towards the east of the network. REMEDHUS consists of four automatically recording meteorological stations and about 20 SM monitoring stations among which more than 10 stations record multi-depth SM. It has already been used for validation of SM content and variability of previous SM products from SMOS and SMAP mission [19,20,66]. The stations are equipped with HydraProbe R (Stevens Water Monitoring Systems Inc., Oregon, USA) and EnviroSMART R (Sentek Pty. Ltd., Stepney, Australia) sensors and provide SM content on an hourly basis at 5 cm, and 25, 50 and 100 cm, respectively, with a nominal accuracy of 0.01 m 3 m −3 . The REMEDHUS network participates in the International Soil Moisture Network (ISMN), an open access database containing global in-situ SM networks that guarantees a standard of datasets for validation and calibration of air-and spaceborne remote sensing products [14]. The homogeneity of the REMEDHUS area in terms of topography and climate provides favourable conditions to compare large-scale satellite imagery with the in-situ measurements. The region is characterized by a semiarid continental climate with an average annual temperature of ∼12 • C. The annual precipitation has been 385 ± 100 mm in the last 10 years and is homogeneous among the network. However, intensity and timing of meteorologic conditions vary over the course of a year. Climate factors in combination with present growing and dormant periods can be used to define SM seasons, which are accompanied by similar patterns of SM recharge and consumption processes. In this work, SM seasons were estimated based on the relationship between PET-the net SM demand if water was potentially available-and precipitation. Results of the comparison of PET and precipitation are given in Section 3.1. Hereby, PET was computed on the basis of the Penman-Monteith equation and was used as a reference to directly estimate the actual crop evapotranspiration [67]. It was obtained from solar irradiance, relative humidity, mean temperature and wind speed. The average of PET over all four meteorological stations was computed. This approach has been shown to accurately estimate evapotranspiration for cropped surfaces in various climates and has been used as the standard method by the United Nations Food and Agriculture Organization.

DTW Technique
DTW is used to study dissimilarities between time series by determining the optimal match between two observations [68]. It is a dynamic programming technique which breaks down a complex problem into an easier solvable sub-problem [48]. Hereby, the initial problem is optimized in a way that the solution of the corresponding sub-problem produces the identical result with a considerable reduction of computational effort. This is particularly beneficial when datasets become large or highly dimensional. The DTW algorithm is based on a local distance measure which contains the distortions and shifts between the analyzed time series. Its main objective is to find the optimal warping path by minimizing a local distance measure between the series within a particular observation period. Two general cases can be pursued-(1) finding a repetitive sub-sequence within a longer sequence and (2) aiming for the optimal global match of two time series within an entire period of observation. Regarding the first case, it was applied in speech recognition and to improve the assessment of electrocardiograms to recognize essential heartbeat patterns within a longer sequence of waveforms [48,[69][70][71]. The algorithm has the advantage to account for a relative warping speed between the time series and considers its local variations. The temporal evolution obtained by DTW can be analyzed to retrieve information about mutual dependencies or causal relationships between the time series. In reality, these time series can represent regularly sampled, continuous geophysical or meteorological parameters which are considered to have a variable time lag over an observed period.

Fundamentals of DTW
The DTW algorithm evaluates two input time series which consist of a sequence X of length M and a second sequence Y of length N with not necessarily the same length as X : The absolute values of two regularly sampled sequences are pairwise compared on the basis of a local distance measure. A criterion to select an adequate distance measure is the dimensionality of the time series [72]. Since only two series are compared, Euclidean distance is used as a local distance measure. The optimal warping path of maximum alignment between two series corresponds to the minimum global cost over the entire period. The global cost gc is defined as the sum of local costs of the L pairwise elements d_l = d(x_m_l, y_n_l) in a local cost matrix along a warping path p: Accordingly, gc indicates the global dissimilarity between time series, that is, the higher the value, the more distinctive the features of the input time series. Its absolute value depends on both the total length and the sampling (temporal resolution) of the time series. Therefore, it is difficult to generally ascribe a physical meaning to this absolute quantity. The warping path can be defined by a set of (X, Y)-tuples that correspond to the L pairwise elements x i ∈ X and y j ∈ Y along its trajectory. It represents the relation between samples of time series X with samples of time series Y: The warping path is implemented by satisfying the following criteria: (i) boundary condition: The boundary condition defines the initial and final states of the warping path. As indicated in (i) the time series are forced to be aligned at the beginning and the end of the observed period. The monotonicity condition ensures causality and particularly refines the path to exclusively advance forward in time. The continuity condition requires the path to go through every sampled point in time. Hereby, a step-size condition is introduced and chosen in a way that both monotonicity is guaranteed and no samples are skipped (continuous path). In a simple version, two consecutive samples either maintain the present time lag (1,1) or change in relative delay or alignment (1,0), (0,1). The optimal warping path is successively obtained by following three steps: Building the accumulated cost matrix C acm (iii) Retrieval of the optimal warping path p opt First, D E is computed from the pairwise local distances d at every point in time: Determination of the optimal warping path from comparison of the global cost of all possible paths is computationally expensive (complexity O(N · M), that is, exponentially increasing cost for an increasing length of time series). Since the optimal warping path is obtained based on a single constant distance matrix D E and satisfies the criteria (i)-(iii) (time series evolve positively monotonous), a successive calculation of the warping path also yields the optimal solution. An accumulated cost matrix C acm is computed from D E by summing up the partial costs given by the corresponding distances of each path section. C acm can be obtained using a simple step-size condition (criterion (iii)) with the elements computed as follows: The left and bottom boundaries in C acm are set to infinity and the first element is initialized with the corresponding value of the bottom-left element of D E . The matrix C acm is then successively calculated from (1,1) to (M,N). An intermediate value at element (m,n) is obtained by summing up the distance in element (m,n) of D E and the minimum accumulated cost of the corresponding element according to the step-size condition. Once C acm is obtained, the optimal warping path is retrieved by backtracking the entries along the 'valley' of the minimum accumulated cost from (M,N) to (1,1). The sum of the distances of all warping path tuples (X, Y) in D E yields the minimum global cost.
Relative changes of global cost between time series at different locations or at different times can be used to describe temporal and spatial variability of the observed parameters. However, being a global quantity over an entire period, it can miss out local information which is needed to understand the behaviour of the warping path and investigate the underlying processes. Therefore, this study is focused on the extraction of a meaningful warping path evolution rather than optimizing the global cost. Geophysical parameters can be retrieved using measurement techniques with differences in data quality as well as temporal and spatial resolution. DTW allows to compare time series of either the same or different parameters which, in addition, can be extracted at different locations. In case the same geophysical parameter is observed and compared at different locations, temporal variations of the warping path may contain spatial variability of the governing processes. Additionally, local warping behaviour between series measured at the same location may reflect their temporal variability. The local information at a particular sample is given by the relative time lag between time series at a specific time. In this study, the time lag evolution between SMOS and in-situ SM observations was investigated. Since the observed parameters have a geophysical meaning, both the evolution and range of the time lag are naturally constrained by their underlying properties and mechanisms. Meaningfulness of the resulting time lag in terms of the reflection of prominent features of the input time series as well as computational efficiency can be improved through customization of the DTW technique (Section 2.2.2).

Customization of DTW
In this section, possible customizations to the DTW technique are explained. They comprise a maximum allowed time lag, the adjustment of the step-size condition and the computation of the dates of onset of pronounced precipitation events. Results of the intermediate steps after application of customized DTW to SMOS and in-situ SM are presented in Section 3.3.

Maximum Allowed Time Lag
The standard implementation of DTW requires the computation of the entire C acm (Equation (5)). With some knowledge about how geophysical observations are temporally related, the range of the warping path can be roughly estimated beforehand. A Sakoe-Chiba band is a global constraint that adjusts a time interval in C acm in which the warping path can freely evolve [49]. It has been widely used in DTW to define the range of maximum allowed time lag (Figure 2a) [68]. After applying a Sakoe-Chiba band, all off-diagonal elements of C acm outside a certain interval around the main diagonal are set to infinity. The evolution of the warping path is limited to a reasonable predefined time lag interval as follows: where lead and delay represent the maximum allowed lead and delay of series Y with respect to series X , respectively. Short-term fluctuations in SM observations typically are responses to certain events (e.g., precipitation). In cases like insufficient temporal sampling, SM fluctuations originated from different events might be incorrectly assigned to the same event. Also, without the Sakoe-Chiba constraint, the warping path can skip important temporal features and, as a consequence, unnaturally long delays are accumulated. SMOS SM is more sensitive to meteorologic conditions. In contrast, in-situ SM observations appeared to be more damped depending on the soil properties [73] and long-term climate factors. Furthermore, recharge of subsurface SM content via infiltration is a gravitationally driven process. Therefore, in-situ SM are assumed to lag behind SMOS observations and hence no lead is expected.

Adjustment of Step-Size Condition
Uncustomized DTW allows the warping path to freely evolve under satisfaction of criteria (i)-(iii). The step-size condition controls the slope of the warping path and defines the sensitivity level to relative changes between observations. It is important to note that the simple version of criterion (iii) corresponds to an uncontrolled slope, which permits the assignment of an infinite number of samples of series Y to one sample of X , and vice versa. An instantaneous accumulation of high time lag between two consecutive samples is unphysical. Conversely, when the slope is controlled, the number of assigned samples per sample is limited and the warping path is forced to advance in time within a maximum defined rate.
The slope control also serves as a low-pass filter to account for highly variable or noisy observations which would lead to overaccumulation of time lag. Additionally, the temporal sampling of the input series affects the time lag evolution for a particular step-size condition. As an example, in case of an uncontrolled slope and fine-sampled series are compared, small fluctuations between the series may already cause unreasonably high variations in time lag. For this reason, it is necessary to control the slope to find a trade-off between the temporal sampling limit, the order of noise level, and the expertise about temporal scale of the processes to be resolved. The 2-3 days revisit time of SMOS defines the sampling limit for application. While near-surface SM time series are generally characterized by fast increase during recharge after precipitation and somewhat exponential decrease during a drying period [39], subsurface SM is mainly governed by intrinsic SM properties and the exposure to continuous atmospheric forcing. Therefore, only pronounced precipitation events lead to notable short-term increase in SM at deeper levels, whereas weaker events cause gradual increase over the long term [74]. The positive and negative slope of the step-size condition can be customized differently. Figure 2b In the step-size condition in Equation (7) the relative rate of reduction and accumulation of time lag between the input series are of 3 and 2 samples per sample, respectively. As an example, step-size interval can be determined by considering the rate of reduction of time lag in case a pronounced precipitation events leads to fast SM increase and the rate of accumulation of time lag during SM decrease with the absence of precipitation.

Determination of Onsets of Pronounced Precipitation Events
Surface processes often affect subsurface SM over long periods of time. However, pronounced precipitation events can lead to a rapid increase in both SMOS SM and in-situ SM down to a certain depth. In this case SM time series are sensitive to the same event. Regarding DTW, the warping path is assumed to be aligned at these particular dates of high precipitation. Subsequently, it further evolves depending on the prominent features as before. A function was developed based on the cumulative sum of precipitation to automatically detect the dates of onset of the k most pronounced precipitation events within a given period of time. To obtain the dates, the precipitation is first discretized with a predefined bin rate with respect to the temporal sampling of the time series. Bin edges are determined by considering an equal amount of precipitation within each bin (fixed bin density) and a variable bin width. Hereby, a small and a wide bin width correspond to a period with a high and a low precipitation rate, respectively. The change from low to high precipitation is indicated by a sharp increase in convexity of the cumulative sum of precipitation. The higher the change in curvature, the more pronounced is a precipitation event. The date of onset of a pronounced event is given at the bin edge, where the differences in adjacent bin widths show a local maximum. Subsequently, the local maxima are sorted to obtain the onsets of the k most pronounced events.
In this work, the number of pronounced events was selected according to the precipitation pattern and the most prominent features of SMOS and in-situ SM along the study period (2016 to 2018). Both SMOS and in-situ SM show significant increase after the occurrence of pronounced precipitation during summer, precipitation events which initiate SM recharge after the mostly dry summer, and heavy precipitation during the wetter winter. The selection of the nine most pronounced events was sufficient to capture the dates of the most prominent features, without overcontrolling the evolution of time lag. The corresponding dates of onset are illustrated in Section 3.3. If applicable, both SM time series are cut at these particular dates (red lines), and DTW can be separately applied to the corresponding sub-sequence, starting each time from alignment of the time series and zero time lag.

Comparison of SMOS and In-Situ SM
Among the different REMEDHUS stations, four stations were selected regarding land use, soil texture and the availability of multi-depth SM measurements (marked in the red boxes in Figure 1 Satellite observations at L-band are known to be sensitive to the SM profile within the first 5 cm profile. In croplands, vegetation root zone is estimated to be within the first 50 cm depth, where most of the root biomass is present [2,75]. For these reasons, in-situ SM measurements at 5, 25 and 50 cm were selected to encompass both near-surface and root-zone SM. Land use, soil type and depth-specific SM-related properties of stations E10, J12, L3 and M5 are shown in Table 1. These properties include field capacity (FC) and wilting point (WP) of the soil bulk at the corresponding depths. FC defines the amount of SM which eventually gets retained subsequent to infiltration, after water has completely drained away (after 2-3 days), and WP estimates the minimum SM content a plant is capable to utilize against the soil matric potential to hold back water in the soil. These parameters were obtained from the sand-silt-clay composition of the soil and vary with soil texture. In case these parameters are vertically heterogeneous, they may indicate differences in SM dynamic range as well as different rates of recharge and dry out. Additional information on the predominant soil type of the uppermost soil horizon at the four study sites were obtained from the Agrarian Technological Institute of Castilla y León (ITACyL) [76].
Since the low spatial variations between SMOS L2 pixels over REMEDHUS found in a previous study were in agreement with in-situ site-specific characteristics [19], the SMOS L4 pixels at 1 km also correspond. However, the level how SMOS L4 SM represents the single stations over the course of a year is not clear. A non-negligible part of SM consumption is due to root-water uptake, which is a particularly important process during the growth of vegetation. Similar to timing and quantifying the amount of precipitation in SM recharge, root-water uptake depends on the land use, which in turn depends on water availability and stage of plant development. Stations M5 and J12 comprise rainfed winter cereals, which is the principal land use of REMEDHUS. Since stations E10 and L3 are located in vineyard, the corresponding SMOS pixels are likely representing a mixed land use. Thus, it is important to know the typical field scale of land use (heterogeneity scale) at the REMEDHUS network, especially at the surrounding of the SM stations. This includes knowledge to which level SM at 1 km resolution can represent point-scale observations and, whether crop-related processes are supposed to be captured in their time series. In this study, the NDVI was used to quantify the variability of land cover in the REMEDHUS network. Temporal variability of land use was studied using 16-day time series composite from MODIS within the three-year study period. High-resolution (10 m) NDVI maps are computed using near-infrared and red bands of Sentinel-2 L1C satellite images with a maximum cloud coverage of 30 % [57]. The following steps were carried out and the results are presented in Section 3.2. The temporal variability of land use was studied using the Terra MODIS NDVI from 2016 to 2018, particularly the 16-day composite at 1 km resolution (MOD13A2 product of collection 6), derived from atmospherically-corrected reflectances at red (band 1: 620-670 nm) and near-infrared (band 2: 841-876 nm). The spatial variability of land use within the corresponding SMOS L4 pixel was investigated using Sentinel-2 NDVI images at 10 m at a representative date within the growing season of winter wheat (21 February 2017), when differences in land use were particularly pronounced. The NDVI was computed from Sentinel-2 L1C reflectances at 10 m in red (band 4: 665 nm) and near-infrared (band 8: 842 nm) with a maximum cloud coverage of 30 %. The mean and variation of NDVI were studied around the stations to understand the level of representativeness of SMOS L4 at 1 km to each of the four in-situ stations.
High-resolution SM from SMOS is compared to in-situ observations at different depths from 2016 to 2018. A three-days average was applied to all SM time series to obtain regular sampling from irregularly sampled SMOS observations. Subsequently, a low-pass filter (Gaussian smoothing σ = 2) was applied to all series to smooth out high magnitude peaks of short-term SM fluctuations and to obtain the comparable trend series. To account for different dynamic ranges of the observations, input time series were calculated by applying min-max normalization to the trend series.

Results
The results of the temporal variability of climate factors in terms of SM seasons are presented in Section 3.1. The estimated levels of spatial representativeness of each station from the analysis of spatial heterogeneity of land cover due to the presence of crop growing seasons are given in  Figure 3 shows the monthly accumulated precipitation recordings and PET estimations, averaged over eight years (2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). PET and precipitation are compared to subdivide a year into three distinctive SM seasons (recharge, utilization and deficit). The recharge season starts as soon as monthly accumulated precipitation exceeds PET around November and coincides with the approximate crop sowing date. While evaporation decreases during the winter months, SM starts to recharge and gets already partly utilized by early plant development and accompanied transpiration. The utilization period starts once monthly accumulated PET exceeds precipitation around mid-February. Water availability from precipitation cannot compensate for losses due to increasing evapotranspiration and vegetation mainly utilizes stored subsurface water. Due to decreasing precipitation and continued water consumption during the growing season, soil starts to dry out. The limit between utilization and deficit is defined to be around mid-June when root-zone SM gets intensively taken up by plants, with interannual variation. Crops are harvested shortly afterwards when water becomes scarcer. Summer months are generally too dry and the associated subsurface SM deficit is not suitable for crop cultivation on a purely rainfed basis. Information about the main features of the SM seasons are summarized in Table 2.   Figure 4 shows the MODIS NDVI time series at 1 km resolution for all stations from 2016 to 2018. Albeit land use differs among the locations of the stations, all NDVI time series generally indicate the phenological characteristics and timing of a winter cereal growing season [55]. After the sowing date around October and subsequent germination, a small green-up peak at the turn of the year represents plant emergence and growth. After a short crop dormancy due to the low winter temperatures, indicated by a local NDVI minimum, rapid plant growth is reinitiated until peak greenness is reached around April. The harvest date is around June/July depending on meteorological conditions and crop maturity. Although total water consumption is high due to the relatively long growing season, summers in the REMEDHUS region are too dry to cultivate summer crops.

Spatial Heterogeneity of Land Coverage
Due to high irradiance and temperature, out of about 500 mm integrated PET over the entire growing season, around 350 mm are estimated to be within the main growth period from May to July (see Figure 3). In this main period the water consumption by root-water uptake is highest. Winter cereal achieves significant root depth in early stage with and an ultimate possible rooting depth above one meter. Thus, the accessibility of deeper soil layers guarantees yield stability even in regions with low rainfall (annual precipitation lower than 600 mm) [77].
To better understand the spatial variability within the SMOS L4 pixel, Sentinel-2 NDVI images at 10 m are shown in Figure 5a-d for the four in-situ stations, respectively. The selected images with full coverage correspond to a representative date within the crop growing period (21 February 2017), when vegetation of the agricultural fields is already sufficiently developed to distinguish a variety of land covers. The typical field scale of cultivation in the REMEDHUS network is significantly smaller than 1 km. Figure 6a-d show the NDVI mean and standard deviation (stDev) of Sentinel-2 NDVI as a function of the square area around the corresponding station. It is worth to note that in all cases, NDVI in a square area of 1 km × 1 km (delta of 500 m) resembles the mixed land use over the entire REMEDHUS network, which was observed to be of 0.32 ± 0.18 at that particular date. Through a visual inspection, the typical NDVI values for vineyard and rainfed cereals were determined to be around 0.15 and 0.45, respectively. Station E10 is located in a vineyard area next to a field with cereal crops. Accordingly, a fast increase of NDVI values for increasing delta and higher variability are given in the close surrounding of E10. Station L3 is more centrally located in vineyard and resembles its typical NDVI of 0.15 up to about 100 m resolution. The increase of mean and variation of NDVI beyond that scale indicates that different land uses are incorporated in the wider surrounding. Station J12 shows an increase in both mean and variation of NDVI for larger resolution, while including crop fields at different stages of development. Regarding both mean and variation of NDVI, among all stations, M5 overall resembles the values at SMOS resolution best and represents similar land use within the entire resolution ranging from a scale of 10 m to 10 km. The high NDVI standard deviation at station M5 at small resolution is due to its location near the border of different land uses. Albeit knowledge of soil texture, spatial variability and distribution of land cover are necessary, additional information of climate factors and the presence of crops are crucial to further understand the plant-moisture relationship and to eventually point out regions with similar SM dynamics. All information about the level of spatial representation of land use of single stations within the corresponding SMOS L4 pixel is summarized in Table 3.  Table 3. Results of the analysis of land use variability in the surrounding of the in-situ locations (delta = 25 m) and at 1 km resolution (delta = 500 m), including an estimated level of spatial representativeness of SMOS and the corresponding in-situ station. NDVI Images (mean values and standard deviation) are evaluated at a representative date during a growing season (21 February 2017).

Interim Results of DTW and Customization
DTW was applied to SMOS and in-situ SM observations and the single steps are explained for station E10 in 2017. For this station, the SMOS and in-situ time series at different depths together with the final results for the entire study period are shown in Section 3.4. The intermediate steps of the application of DTW to SMOS and in-situ 5 cm time series are illustrated in Figure 7 as an example. This includes the computation of the local distance matrix D E from the normalized trends, the accumulated cost matrix (C acm ) with a particular step size ∈ [−∞,3] and a maximum allowed time lag of 10 samples (30 days). The time lag evolves along the optimal warping path for the given customization, indicated by a red line. Figure 8 illustrates the customization of DTW for station E10 at 50 cm depth in 2017. The particularly dry year with intermittent pronounced precipitation events in summer clearly led to short-term recharge at 50 cm depth. Figure 8a shows the corresponding normalized trends in combination with the three-days accumulated precipitation. SMOS SM increases at almost all precipitation events, while in-situ SM is only sensitive to the most pronounced events, such as to those in the beginning of the July, September and November, when the soil was initially dry. It can be seen that SMOS SM decreases several days after the absence of precipitation, whereas in-situ SM is only characterized by gradual drying. Figure 8b shows the time lag evolution, corresponding to a maximum allowed time lag of 60 days under customization of the negative slope and fixed maximum positive slope of two samples per sample, respectively. Irrespective of the customization of the negative slope, the time lag runs into the maximum allowed value between May and mid-July. Thus, the time lag does not reflect the observed alignment of the input series at the onsets of the most pronounced precipitation events in the beginning of May and July, when both trend series resemble clearly. For a finite negative slope of −2 or −3, the time lag gets further trapped at high accumulated values until mid-October. The customization with an infinite negative slope performs best, but is not sensitive enough to lead to a prompt reduction of time lag, if necessary. To enhance the sensitivity to common features of the input series, the alignment of time series are forced at the onset dates of the most pronounced precipitation events. This results in zero time lag at that particular dates and allows the warping path to evolve naturally after alignment is forced. Consequently, the four pronounced events in 2017 were selected when SMOS and in-situ trend series show clear alignment (indicated as red lines in Figure 9 in 2017). DTW was subsequently applied to the resulting sub-periods in 2017, which are limited by the four onset dates. The corresponding time lags of the sub-periods were eventually merged. Figure 8c shows the time lag evolution after forced alignment, including different customization of the positive slope and infinite negative slope. The dates when alignment is forced are marked by red crosses. It can be seen that the main features of the input series are reflected in the evolution of time lag. Features which are only visible in SMOS SM (e.g., SM increase after a rain event in mid-October) do correctly lead to further accumulation of time lag. Therefore, the automatized implementation to force the alignment increases the sensitivity of DTW; the customization with an infinite negative and finite positive slope showed the overall best performance for all depths.

Final Results of the Evolution of Time Lag
In this section, the results of time lag evolution between SMOS and in-situ SM time series at the depths of 5, 25 and 50 cm of the representative stations M5, E10 and J12 from 2016-2018 are presented. In all cases, the time lag is determined from customized DTW using a step size ∈ [−∞,2] and forced alignment of the time series at the onsets of the nine most pronounced precipitation events within the study period. Customization of a maximum allowed time lag was not applied. The results for each station are given in Figures 10-12. The time lag evolution of the corresponding time series are displayed below. The onset dates at which time lag gets reseted to zero are depicted by red lines in the distribution of precipitation along the study period (Figure 9, previous section), and by red crosses at the resulting time lags (Figures 10b, 11b and 12b). The results were assessed in terms of SM time series, the previously defined SM seasons ( Table 2 in Section 3.1), and the level of spatial representativeness of the SMOS observations ( Table 3 in Section 3.2). The best results were obtained for stations M5 and E10, where land use between SMOS and in-situ station is well-represented. For these particular stations, the following general observations can be made regarding amount and timing of accumulation of time lag.  Figure 9), when dry conditions are followed by abundant rains. This leads to an increase in subsurface SM and to more frequent alignment of the time series within the recharge season. Thus, almost no time lag is observed in the period from October to January and both subsurface and SMOS SM follow the same pattern. A small time lag implies that SM changes mainly depend on shortly preceding atmospheric conditions, which is in agreement with the fact that in this case the increase of SM in the recharge process cannot be referred to a low initial SM content in the beginning of the recharge season. After sufficient recharge occurred at all depths, continuous precipitation and moderate consumption are barely changing the SM content. As a consequence, SM at root zone becomes widely insensitive to changes at the surface, which results in less accumulation time lag. With increasing water consumption and evapotranspiration during utilization, SM gradually starts to dry. This dry out gets interrupted by intermittent precipitation, which results again in low accumulation of time lag. In the beginning of the deficit season, SM decreases substantially due to both absence of precipitation and strong root-water uptake. After crops are harvested, less-intense, but continuous SM decrease is recorded close to the WP, where water becomes unavailable for roots. From now on, evapotranspiration becomes too high for effective infiltration and less frequent rain events marginally influence SM. These circumstances lead to high accumulation of time lag up to the length of the entire deficit season, especially for deeper observations, and the actual SM content is more influenced by initial SM conditions. 2017 was an anormal year regarding both timing and amount of precipitation with overall dry conditions (annual rainfall only half the usual average), which had consequences for the time lag evolution. Firstly, due to the absence of precipitation, high accumulation of time lag could be already observed in March. Furthermore, the two occasional strong rainfall events in the deficit season led to SM recharge and hence time lag was reduced.
Regarding soil texture, clay-rich soils with high FC are characterized by less pronounced SM changes in comparison to sandy soils. This can be seen at station J12 consisting of clay-rich soils, where SM recharge at 50 cm only occurs at the very end of the recharge season. It is important to note that the presence of crops can influence the time lag significantly as well. This can be especially seen at station J12 towards the end of the main growing season. In the absence of precipitation around July 2016 and 2018, strong root-water uptake leads to a fast decrease of root-zone SM at 25 and 50 cm depth, which is unusual for clay-rich soils. SM at 5 cm is less utilized and remains more sensitive to atmospheric forcing. This has a reverse effect on the time lag showing the highest accumulation for 5 cm SM along the depth profile. Station L3 was categorized as the least representative among the four stations. Due to limited validity, results of this station are presented, but the time lag evolution is not further analyzed in detail. The time lag within the estimated SM seasons are summarized in Table 4 for each station and depth. Mean and maximum values in each period are averaged over the representative years 2016 and 2018.

Discussion
In this study, the transient SM dynamics at near surface and root zone was investigated by quantifying the typical response time between SMOS satellite and in-situ observations. A DTW approach was customized to obtain the time lag between their trend series-which is assumed to be a dynamic quantity rather than a single constant-describing the relationship of their common prominent features. The aim was to assess whether repetitive patterns of alignment and accumulation of time lag can be attributed to SM-related processes on the basis of additional information about climatic factors, soil properties and heterogeneity of land use among the study area. Differences in the evolution of time lag could be referred to specific SM seasons, and variations due to the level of spatial representativeness of land use between SMOS and in-situ observations. The results suggest that causality information about processes can be reflected in the time lag. In former studies, the response time of subsurface SM observations to satellite time series were assumed to increase along the depth profile [18,38,39,41]. However, during the main growing season of cereals, root-water uptake led to an inversion of the typical response time, showing longer time lags at shallow depths. Further investigation is required to quantify the observed phenomena. In this work, a semiarid region was investigated, consisting mainly of a seasonal continental climate with cultivated winter cereals. The scope of this work was not to quantify causality of targeted variables. To do so, the evolution of the time lag has to be evaluated regarding different climates, topographic characteristics and land covers to understand how dependencies and certain processes are eventually reflected in its local behaviour.
SM dynamics depend on atmospheric conditions, local soil properties, and the initial SM state. The availability of precipitation data is crucial to estimate both amount and timing of SM infiltration and to distinguish between rain events of different intensities. In case subsurface SM reaches values close to saturation or dry-end of the soil, the contribution of atmospheric conditions is limited and the time lag becomes less representative. During the deficit season for example, the SM release towards the atmosphere through evaporation is particularly high. Consequently, SM at near surface increases in response to moderate-intense rainfall events, but water is not infiltrated into deeper layers and hence, SM dynamics at 50 cm are not captured in the satellite observations. After intensive recharge, water can be retained in the soil for weeks and subsurface SM values are reaching their maximum close to saturation. Thus, even in cases of ongoing rains, no significant changes of subsurface SM are expected, and the time lag becomes insensitive to the recorded SM changes at the surface. The matching process between the time series was improved by separating the study period into sub-sequences by determining the dates of pronounced precipitation events, to which both time series are sensitive to. Conversely, wavelet transforms of the time series can be analyzed to relate the prominent features of the observations and to approximate representative sub-sequences. This approach could be used as an alternative to customization prior to DTW to improve the matching process by considering the local structural information of the observations. The time lag evolution can be used to assess whether remotely-sensed SM is in agreement with the corresponding ground-truth acquisitions. Studies have shown that satellite SM showed a clear bias in both mean value and amplitude of fluctuations in comparison to ground-truth SM [20,66,78]. In-situ observations generally show a smaller dynamic range than satellite observations, due to limitation of intrinsic SM-related soil properties [79]. The more a pixel is spatially extended, the more constrained is the dynamic range of the area-averaged satellite observation. Biased values can have non-negligible seasonal variation [80]. Furthermore, penetration depth of satellite measurements is not constant, but depends on the attenuation of the signal due to changes in SM content in the uppermost soil layer and hence affects its validation [81]. As a result, the choice of an appropriate metric should be based on the nature of the retrieved variables (dynamic range and fluctuations) including advantages and drawbacks for a particular application [82]. The time lag evolution provided in this paper relies on relative behaviour of the input series rather than on absolute values only and has the advantage of accounting for temporal shifts and non-stationary biases between time series.
In former studies, estimates of root-zone SM by applying an exponential filter to the remotely-sensed observations under the optimization of a depth-specific characteristic time length over a period showed good results. However, in critical periods of sudden SM changes, the assumption of a constant time length may lead to oversmoothing in the estimated values depending on the weighting criteria of preceding SMOS measurements. Knowledge of the climate factors, the presence of growing seasons and the accompanied root-water uptake as well as the level of spatial representativeness of the satellite observation to field scale can be added to infer subsurface SM more accurately. An upfront determined time lag, which is tailored to a corresponding target region and period, can be utilized as a measure similar to the characteristic time length to improve performance in the inference of subsurface SM from satellite observations.

Conclusions
This study presented a methodology to estimate the response time between SMOS satellite and point-scale in-situ SM observations at near surface and root zone. The response time is assumed to be highly dynamic, showing seasonal variability and variations among different soil properties and land uses. The DTW technique is supposed to cope with non-linear effects in the time series to provide a time lag which captures the relative SM changes between the observations. DTW is capable of relating observations based on the relative behaviour of their common prominent features and was used to quantify the response time. SM observations at 5, 25 and 50 cm of four representative stations within the semiarid REMEDHUS SM network were selected regarding their land use and soil properties.
The focus was on investigating the prevailing SM-related processes and local variabilities to understand their impact on the resulting time lag. Climate factors are temporally variable and show spatial variations only at a relatively large scale, whereas the presence of crops and the related root-water uptake during growing seasons occur at smaller field scale, causing SM to be spatially and temporally variable. Based on the comparison of long-term precipitation recordings and PET estimations over REMEDHUS, the classification of periods (recharge, utilization and deficit) with similar patterns and the distinction of the prevailing processes among these three seasons support the interpretation of time lag. MODIS NDVI time series at 1 km revealed the phenological cycle of the predominant land use of the study area, which mainly consists of rainfed cereals. Thus, the spatial heterogeneity of land cover was analyzed from Sentinel-2 NDVI images at 10 m to estimate the level of spatial representativeness of the in-situ station within the corresponding SMOS L4 pixel at 1 km during the main growing season. The spatial scale of cultivated fields is typically much smaller than 1 km, which led to non-negligible differences in heterogeneity among in-situ stations with respect to the SMOS L4 pixels.
DTW was applied to the normalized trend series of SMOS and in-situ observations from 2016 to 2018 and the evolution of time lag was obtained. Customization was applied to improve the accuracy of the technique to ensure that common prominent features of both input series are reflected in the time lag. Positive and negative slopes, which correspond to the rate of accumulation and reduction of time lag, were controlled to consider the natural behaviour of drying periods and wetting conditions. In the case of the anormal year 2017, intermittent heavy rains within drier conditions in summer incorrectly led to an unreasonable accumulation of time lag, whereas time series showed clear alignment. To account for this phenomenon, a method based on the curvature of the cumulative sum of precipitation was implemented to determine the onsets of pronounced precipitation events. After these particular events, when alignment of the time series was forced and the time lag was reseted to zero, the time lag continued to evolve naturally.
The following conclusions can be drawn regarding the evolution of time lag at the particular stations. Finer-textured soils result in higher time lag. The deficit season clearly shows the highest time lag. When abundant rains prevent long-term dry out during recharge and utilization seasons, essentially less time lag is accumulated. In most cases, time lag increases with increasing depth. This is given for the two most representative stations M5 and E10. An inverse behaviour could be observed at station J12, where strong root-water uptake was present during the main growing season. Hereby, SM at root zone was consumed rapidly, whereas soil at shallow depths dried out more slowly, which resulted in less accumulation of time lag.
This study focused on the link between SMOS and in-situ SM series to characterize the involved SM-related processes by observing a total of four representative stations. In a next step, the variability of the time lag can be further investigated by comparing SM stations with similar climate factors, to better point out in which cases integrated climate parameters and information about land cover are necessary. Once the behaviour of time lag is better understood, general conclusion can be drawn if the approach is also applied to different SM networks located in different climates. Funding: The project that gave rise to these results received the support of a fellowship from "la Caixa" Foundation (ID 100010434). The fellowship code is LCF/BQ/DI18/11660050. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 713673. This study was also funded through the award "Unidad de Excelencia María de Maeztu" MDM-2016-0600 and by the Spanish Ministry of Science and Innovation through the projects ESP2017-89463-C3-1-R, ESP2017-89463-C3-2-R and ESP2017-89463-C3-3-R.