Field-scale apparent hydraulic parameterisation obtained from TDR time series and inverse modelling

Due to the large heterogeneity in the hydraulic properties of natural soils, estimation of field-scale hydraulic parameters is difficult. Past research revealed that data from accurate but small-scale laboratory measurements could hardly ever be transferred to the field scale. In this study, we explore an alternative approach where apparent hydraulic properties of a layered soil profile are directly estimated from hydraulic inverse modelling of a time series of in situ measured soil water contents obtained from time domain reflectometry. The data covered a one-year period with both wet and dry soil conditions. For the time period used for inversion, the model is able to reproduce the general evolution of water content in the different soil layers reasonably well. However, distinct drying and wetting events could not be reproduced in detail which we explain by the complicated natural processes that are not fully represented in the rather simple model. The study emphasises the importance of a correct average representation of the soil-atmosphere interaction.


Introduction
Soil hydraulic and transport properties are key parameters required for modelling water and solute movement in soils.The estimation of their effective values at the field scale remains a challenge for hydrology, however.Original concepts envisaged soil as an essentially uniform medium.In this context, properties could be determined by extracting small samples and by performing appropriate measurements on them.Initially, these aimed at the direct determination of the properties of interest (e.g., Topp and Miller, 1966;Klute, 1986).It turned out, however, that these methods yield rather inac-Correspondence to: U. Wollschläger (ute.wollschlaeger@iup.uniheidelberg.de)curate results (Flühler et al., 1976).This led to the development of inverse methods where a given parameterisation is adjusted such that modelled results are in optimal agreement with the corresponding measurement.Comprehensive reviews are provided, e.g., by Hopmans and Šimůnek (1999) and Vrugt et al. (2008).Significant effort has been spent to improve the methods from the original One-Step Outflow (Parker et al., 1985), through Multi-Step Outflow (van Dam et al., 1994), all the way to evaporation experiments (e.g., Šimůnek et al., 1998;Schneider et al., 2006).As a result, the hydraulic properties of a soil sample can now be determined rather accurately and over a wide range of hydraulic states.Compared to field measurements, these methods are very time consuming, however, hence expensive.
A more severe problem, whose implications were not realised immediately, arose with the insight that the subsurface is heterogeneous (Nielsen et al., 1973) and even worse that this heterogeneity is of a hierarchical nature (Gelhar, 1986;Cushman, 1990;Vogel and Roth, 2003).Indeed, a number of studies reported disagreement between hydraulic parameters determined by laboratory methods and parameters estimated with direct or inverse methods using field data (e.g., Dane and Hruska, 1983;Ritter et al., 2003;Wöhling et al., 2008).With this, the question of effective material properties emanated, an issue that has been elucidated to great depths for the linear dynamics of groundwater flow and transport (e.g., Dagan, 1986;Gelhar, 1986;Neuman and Di Friderico, 2003) but whose discussion has been rather anecdotal for the highly non-linear regimes of soils.The general concept behind "effective models" is to find a uniform (macroscopic) representation of a heterogeneous (microscopic) reality such that the macroscopic dynamics mimics the averaged microscopic dynamics.This typically entails the deduction of an appropriate differential equation, the model, and of appropriate parameters.An example is Richards' equation, together with soil water characteristic and hydraulic conductivity function, as an effective macroscopic representation of Stokes flow at the Published by Copernicus Publications on behalf of the European Geosciences Union.
U. Wollschläger et al.: Field-scale parameterisation from TDR time series pore scale.Model and parameters are called "effective" if they have been shown to be faithful macroscopic representations.They are called "apparent" if this is merely postulated, as is the case in this article.
For the time being, until solid theoretical solutions are available, the inverse estimation of hydraulic parameters from in situ measured state variables is a pragmatic path to follow.Examples of readily available state variables are volumetric water content and matric potential which can be measured continuously by time domain reflectometry (TDR) and tensiometry, respectively.The advantage of using in situ data is that the sensors experience all processes affecting the measured state variables in their natural environment, implicitly including interactions between different soil layers and across scales.This can hardly ever be achieved in laboratory measurements.
A number of inverse estimations of hydraulic parameters in field soils have been conducted using data from infiltration and drainage experiments in lysimeters (Dane and Hruska, 1983;Abbaspour et al., 1999;Sonnleitner et al., 2003;Mertens et al., 2006).In addition to state variables, these experiments also yield data on leachate volume which further constrain the model.
So far, only a few studies aimed at the inverse estimation of hydraulic parameters from state variables measured directly in natural soil profiles (Vereecken et al., 2008).Lehmann and Ackerer (1997) used pressure head data measured in one single soil layer together with precipitation and evapotranspiration measurements from a one year period to determine hydraulic parameters of a two-layer soil profile.Abbaspour et al. (2000) ran an irrigation and drainage experiment in a layered field soil and estimated the hydraulic parameters of the profile from transient data of water content and matric potential.Doing measurements in a rather complicated soil profile they found that a fairly correct representation of soil structure in the model is essential to infer a reasonable set of parameters from the inverse simulations.Jacques et al. (2002) estimated hydraulic and solute transport parameters of a layered field soil under natural rainfall conditions using time series of water content, electrical conductivitiy and pressure head.In their experiment, evapotranspiration was eliminated by a thin gravel layer placed on the bare soil surface.Ritter et al. (2003) conducted an inverse estimation of hydraulic parameters from a time series of measured water contents from different depths of a temporarily irrigated, layered field soil on Tenerife.Their simulations suffered from ill-posedness of the inverse problem however which did not lead to a global solution.Ritter et al. (2003) attributed this to the large number of hydraulic parameters to be estimated and suggested a need for more experimental data like pressure head and/or outflow data to additionally constrain the inversion.Recently, Wöhling et al. (2008) compared three different multiobjective optimisation algorithms for the estimation of hydraulic properties from pressure head data measured in a layered field soil on the Spydia site, New Zealand.Com-pared with forward simulations using hydraulic parameters derived from laboratory experiments all inverse models were more successful in estimating parameter sets that reproduced the measured field data.Differences in some of the parameters derived from the three optimisation algorithms were explained by the occurence of local optimal solutions and the high dimensionality of the inverse problem combined with the applied multiple objectives.
Most of the previous studies either operated in a rather narrow range of hydraulic states with approximately constant water contents, hence in a quasi-linear (meaning that the dynamics can be linearised) regime, and under fairly wet conditions.Even though they were conducted in natural soils, boundary conditions were often controlled in some way, either by eliminating evapotranspiration or by applying artificial irrigation.
The objective of our study is to estimate the apparent hydraulic parameters of a layered soil profile under completely natural and highly variable forcing.To this end we explore the numerical inversion of a time series of in situ measured water contents over a time period of one year.The water content in the topsoil shows a rather intense temporal dynamics which is expected to stabilise the inversion on the one hand but also tends to invoke processes that are not incorporated sufficiently well in most models like, e.g., evolution of cracks under dry conditions, preferential flow or an accurate representation of evapotranspiration.

Site description and instrumentation
Experiments were conducted at the Grenzhof Test Site (49 • 25 N, 8 • 37 E) near Heidelberg, SW-Germany.The site is a former agricultural field which, since the beginning of the experiments in spring 2003, is covered with grass that is regularly cut to a height of a few centimeters.According to the USDA-Soil Taxonomy (Soil Survey Division Staff, 1993) the overall soil texture can be classified as sandy loam (Fig. 1).The topmost 0.28 m consist of a humous plough horizon that is underlain by a rather homogeneous sandy loam.At 0.82 m, we found a sharp transition to a more dense sandy loam including aggregates of secondary minerals.Further down, the clay content of this horizon increases continuously.Below 1.44 m, there follows a transition to a thin layer of gravels within a loamy matrix.It is underlain by fluvial gravels embedded in a sandy matrix starting between 1.54 m and 1.65 m depth.
The site is equipped with an automatic weather station that continuously measures precipitation, air temperature and relative humidity, barometric pressure, wind speed and direction, and incoming and outgoing solar and far infrared radiation (Table 1).Precipitation, air temperature and relative humidity are measured at 2 m height.The radiation sensor Hydrol.Earth Syst.Sci., 13,[1953][1954][1955][1956][1957][1958][1959][1960][1961][1962][1963][1964][1965][1966]2009 www.hydrol-earth-syst-sci.net/13/1953/2009/ is mounted at a height of 2.82 m.Wind speed is measured at 3.18 m above the ground surface.For calculation of reference evapotranspiration, wind speed data were converted to values for 2 m height using the logarithmic wind speed profile formula given by Allen et al. (1998, p. 56).In addition, volumetric soil water content and soil temperature are recorded in three nearby profiles.Meteorological parameters and soil temperature are measured in 10-min intervals.TDR measurements of soil volumetric water content are recorded hourly.
Soil temperature probes are installed in a soil profile at depths of 0.12 m, 0.17 m, 0.265 m, 0.47 m, 0.67 m, 0.865 m and 1.64 m.Volumetric soil water content is measured with TDR using standard 3-rod probes (CS610, Campbell Scientific, Logan, UT; rod length: 0.3 m, rod diameter: 0.48 cm, rod spacing: 2.2 cm).They are installed in two adjacent soil profiles in 0.13 m, 0.72 m and 1.41 m depth (first profile) and 0.30 m, 0.63 m, 0.92 m and 1.16 m depth (second profile).Before installation, TDR probes were calibrated for water content estimation with measurements in water and air.Bulk dielectric permittivities ε c [-] were derived from the measured L a /L values using where L a [m] is the apparent length of the TDR probe rods which varies with water content, and L [m] is the real rod length (Campbell Scientific Inc., 2004).Volumetric water contents θ [-] were calculated using the CRIM (Complex Refractive Index Model) formula , and ε a [-] are the dielectric permittivities of the soild matrix, water and air, respectively, and φ [-] is porosity.The dielectric permittivity of water (ε w ) was temperature corrected according to Roth et al. (1990), and ε s was set to an estimated constant value of 5.The porosity of each soil layer was determined gravimetrically from the dry weight of undisturbed volumetric soil samples taken during a profile excavation.

Hydraulic model
Our model is a highly simplified representation of reality and neglects processes like preferential flow, hysteresis The investigation of the sensitivity of the inversion to all these factors is beyond the scope of this article.In our simple approach, we assume in the following that the model setup using Richards equation, the Mualem-van Genuchten parameterisation, and the chosen material model and lower boundary condition is applicable to describe the general water movement in the soil profile.Clearly, this is a design decision and there may exist various other models which may describe the observed measurements equally well or even better.
Simulations were run using the software package of the numerical model HYDRUS-1D, Version 2.0 ( Šimůnek et al., 1998b) in the forward mode accompanied by external routines for inversion which were programmed in MATLAB (The MathWorks, Inc., Natick, MA).

Governing equations
Slow vertical water movement in a one-dimensional, unsaturated rigid porous medium is described by the Richards equation (Jury et al., 1991) where θ is the volumetric water content The unsaturated hydraulic properties can be parameterised by the Mualem-van Genuchten model (Mualem, 1976;van Genuchten, 1980): where θ r [-] and θ s [-] are the residual and saturated water contents, respectively, α [m −1 ] and n [-] are shape parameters, K s [m d −1 ] is the hydraulic conductivity at saturation, and [-] is the effective saturation.
Root water uptake can be described by the model of Feddes et al. (1978) where the sink term S is defined as with β(h) [-], a stress response function for root water uptake (0≤β≤1) and S p [d −1 ] the potential water uptake rate ( Šimůnek et al., 1998b).

Material model
The material model was set up employing the textural information of the soil profile shown in Fig. 1.The partitioning into different layers was adapted to the observed boundaries of soil horizons (Sect.2.1) and the availability of TDR probes installed in the profile.We chose five uniform layers with a constant vertical discretisation of 0.01 m.The first layer represents the humous plough horizon.It is followed by the homogeneous sandy loam located between 0.28 m and 0.82 m.Soil layer 3 is modelled as two separate layers in order to account for the continuous increase in clay content with depth.
In fact, this textural transition implies also a continuous alteration of the hydraulic properties which is not represented in our material model, however.For the deeper section of the profile we simpified the model by adding the small layer of gravels with loamy matrix that occurs below 1.44 m depth to the fourth model layer.Additionally, we inserted a further layer at the lower end of the profile (1.55 m to 4.00 m) for representing the coarse grained gravels and sands below.This transition is expected to act as a capillary barrier which is assumed to be hydraulically relevant.Due to the high fraction of gravel and stones no further TDR probes could be installed below 1.44 m in the soil profile.
According to the information from the texture analysis we used standard parameters for sandy loam (Carsel and Parrish, 1988) provided in the HYDRUS-1D, software package Version 3.0 ( Šimůnek et al., 2005) database as a first reasonable initial estimate for the Mualem-van Genuchten parameters of the upper four model layers.For each parameter, bounds were used to limit the inversion to a reasonable parameter range (Table 2).To investigate the uniqueness of the estimated parameters, inverse simulations were repeated using loam, silt loam and loam (Carsel and Parrish, 1988) as initial parameter guesses.In all runs, the standard values for the water content at saturation θ s of each layer were replaced by the measured porosities which were assumed to be equal to θ s .For the coarse grained sediments of model layer 5 we used hydraulic parameters of sand given by Carsel and Parrish (1988).Since no TDR measurements in this layer are available these initial values were kept constant during all simulation runs.

Initial and boundary conditions
Estimating the initial pressure head profile for a layered field soil under non-steady state conditions is difficult since the exact distribution of the soil moisture content over depth is not known.Hence, for a first guess, pressure heads for day 777 (the first day in the simulations, Fig. 2) were calculated from the inverse of Eq. 5, h(θ), for the water contents measured at the different TDR probes.These point values were then interpolated linearly for each model cell between the installed sensors.We chose to interpolate pressure heads instead of the measured water contents since water contents Table 2. Example set for initial hydraulic parameters (here sandy loam for layer 1...4 and sand for layer 5 (Carsel and Parrish, 1988)) and parameter bounds used in model layers 1...4 for the inversion.Mualem-van Genuchten parameters of model layer 5 were kept constant throughout the inversion.are discontinuous across layer boundaries and should not be interpolated under such conditions.The pressure heads calculated for the upper-and lowermost TDR positions were extrapolated to the layer boundaries of model layer 1 and 4, respectively (Fig. 3).For the pressure head condition in model layer 5 we assumed hydrostatic equilibrium conditions for initialisation of the model.This, however, leads to an unphysical jump at the layer boundary.To allow for the relaxation of the initial pressure head profile we used a pressure head [m] Fig. 3. Initial pressure head condition (black line) and equilibrated pressure head profile at the end of the spin-up phase (red line) examplified for sandy loam considered as initial parameter set for the inverse simulations.Black dots indicate positions of TDR probes whose water contents were used for interpolation of the pressure head profile.The pronounced peak near the surface results from root water uptake.
spin-up time of 45 days and started the inverse parameter estimation afterwards.The pressure head profile at the end of the spin-up phase is displayed in Fig. 3 as well.
In order to account for potential upward and downward fluxes caused by the capillary barrier at the transition between the fine grained sediments of model layer 4 and the coarse grained gravels and sands of model layer 5, the lower boundary of the model was realised by a groundwater table (h=0 m) at 4 m depth.With the simulated water table some −2.5 m/α≈36 scale heights away from the soil layer above and the saturated hydraulic conductivity of the fifth model layer orders of magnitude higher than the mean water flux, the influence of the artificially set water table is deemed negligible for all practical purposes.In reality, the water table is some 10 m below ground surface.For the transient upper boundary condition we used daily sums of precipitation and daily average values of reference evapotranspiration (Fig. 2).Reference evapotranspiration j ET 0 [mm d −1 ] as part of the upper boundary condition was calculated from daily average values of air temperature, relative humidity, wind speed and net radiation using the FAO Penman-Monteith equation (Allen et al., 1998, p. 24) where R n is the net radiation at the crop surface T is the air temperature [ is the slope of the vapour pressure curve [kPa • C −1 ], and γ is the psychrometric constant [kPa • C −1 ].Since, on a daily basis, the soil heat flux G can be assumed to be relatively small compared to R n (Allen et al., 1998) it was neglected in our calculations of j ET 0 .The FAO Penman-Monteith equation yields a reference evapotranspiration assuming a "hypothetical crop with an assumed height of 0.12 m, with a surface resistance of 70 s m −1 and an albedo of 0.23, closely resembling the evaporation from an extensive surface of green grass of uniform height, actively growing and adequately watered" (Allen et al., 1998, p. iii).
The FAO (Allen et al., 1998, chapters 6 and 7) recommends to account for the difference in evapotranspiration between the reference grass and the site specific crop by invoking a crop factor κ [-].Two different approaches for estimating the crop factor are offered.The first one uses a single crop coefficient κ which scales the reference evapotranspiration such that j ET c =κj ET 0 .The second approach accounts for crop transpiration and soil evaporation separately by splitting κ into a basal crop coefficient κ cb and an evaporation component κ e : κ=κ cb +κ e .The single crop coefficient represents kind of a time-averaged reference evapotranspiration and neglects short-term daily fluctuations while the dual coefficient approach is more sensitive to single wetting events.However, the dual coefficient approach requires information about variables which depend on soil texture such as estimates for field capacity, wilting point and influence depth of evaporation which may vary over considerable ranges and which depend on the hydraulic soil properties which we whish to estimate during the inversion.For both approaches, tabulated values for different crops and growing stages are provided by Allen et al. (1998) but it is still recommended to adjust the coefficients to site-specific conditions.The application of the crop factor still assumes standard conditions meaning that there are no limitations in soil water availability.Hence, influences of soil water shortage, which are accounted for separately in the root water uptake function of HYDRUS-1D, are not considered when using j ET c .
Since the Grenzhof Test Site is covered with grass, in this study, reference evapotranspiration j ET 0 was initially assumed to be equal to potential evapotranspiration.Following Campbell and Norman (1998, pp. 249-250) potential evaporation j E 0 can be calculated using j E 0 =τj ET 0 , and potential transpiration j T 0 is derived from j T 0 =(1−τ )j ET 0 .The fraction of incident radiation τ [-] which reaches the soil surface and is not intercepted by the canopy is derived from where k is the extinction coefficient [-] which was set to a constant value of 0.398 (Ritchie, 1972), and L t [-] is the leaf area index.L t was calculated according to Menzel (1997, p. 34) who determined an empirical formulation for the leaf area index of grassland where h g [m] is the height of the grass.Since we use the FAO Penman-Monteith reference evapotranspiration j ET 0 (Eq.9), the height of the grass was set to a constant value of 0.12 m for calculating L t , even if h g at the Grenzhof Test Site was lower on average.
Root water uptake was simulated employing the model of Feddes et al. (1978).We used root water uptake parameters for grass according to Taylor and Ashcroft (1972), both as provided in HYDRUS-1D.For the rooting depth only a rough estimation from visual inspection of a soil profile was available.Grass roots were densely distributed within the top 0.08 m and decreasing rather sharply below.In the model this was realized by applying a homogeneous root distribution over the first eight centimeters and then decreasing it by a factor of two for each consecutive model cell of 0.01 m thickness over the following 0.06 m.The impact of this parameter was explored by running simulations with various homogeneous rooting depths for the upper dense section and the same decrease in root distribution below.The absolute value of the minimum allowed pressure head at the soil surface required by HYDRUS-1D was set to 150 m.

Inverse parameter estimation
Inverse estimation of the hydraulic parameters was conducted by choosing the objective function where θ meas and θ sim are measured and simulated volumetric water contents, respectively, N is the number of observations, and b is the parameter vector.For the inversion we applied a trust-region solver.In contrast to standard optimisation algorithms such as Levenberg-Marquardt (Marquardt, 1963), trust-region methods do not operate on a single point in parameter space, but approximate the local neighbourhood, and minimise over this region.This increases the stability of the algorithm and makes it less susceptible of getting  r , s ,  , n , K s ; ; h z

HYDRUS1D forward simulation update
Fig. 4. Workflow of the inverse modelling.The initial parameters p i0 for model layers 1...4 were taken from the HYDRUS-1D database (Carsel and Parrish, 1988).In addition to the Mualemvan Genuchten parameters θ r , α, n and K s the initial pressure head profile h(z) was updated subsequently for every single inversion step.Simulations were conducted with and without inverting the crop factor κ.
trapped in local minima.The algorithm we used is described in Coleman and Li (1996).An implementation is provided by the MATLAB optimisation toolbox, which uses preconditioned conjugate gradients (Hestenes and Stiefel, 1952) to solve the linear equation system for each step.
An overview of the workflow of the inversion procedure is shown in Fig. 4. Inverse estimation of hydraulic parameters for model layers 1 to 4 was done by optimisation of the parameter vector b={p 1 , ..., p 4 ; κ} with p={θ r , α, n, K s }.In choosing these four parameters instead of the five implied by Eqs.(5-7), we account for the fact that hydraulic experiments only depend on θ=θ s −θ r , not on θ r and θ s individually.Operationally, θ s is arbitrarily set to the measured porosity, θ r is eastimated in the inversion.In total, this leads to a 17-dimensional optimisation problem.
The inverse simulations were done by using HYDRUS-1D in the forward mode as a module of the parameter optimisation which was implemented in MATLAB.The first inversion run was conducted using the inital guesses for the hydraulic parameters of each layer and the previously calculated initial pressure head profile.After each inversion step, the optimised Mualem-van Genuchten parameters were updated in the forward model.Additionally, the intitial pressure head profile was recalculated and updated in the forward model using the measured water contents from the different TDR probes together with the optimised parameters from the actual inversion step.This way, we aim to improve the initial condition to approach the real initial pressure head profile and furthermore facilitate the spin-up of the model.
For minimisation of the objective function we used daily average values of measured water contents from one of the soil profiles (see Sect. 2.1) of the TDR probes installed at 0.63 m, 0.92 m and 1.16 m depth (Fig. 2).Data from the TDR probe at 0.30 m was neglected since this probe is located very close to the first layer boundary and measured water contents are expected to be a mixture of water contents in the first and in the second layer.For the estimation of hydraulic parameters of the uppermost layer we used water contents measured in the adjacent soil profile at 0.13 m depth assuming textural homogeneity over this short distance of 2.8 m.Recent ground-penetrating radar ground wave measurements from the Grenzhof Test Site (A. Lodde, unpublished data, 2009) show almost constant near-surface soil water contents along distances of tens of meters.Hence, we presume that using the water contents from the nearby profile for inverting the hydraulic parameters of the uppermost model layer is justified.
We used measurements from the time interval ranging from 15 Febuary 2006 to end of January 2007 (Fig. 2).The data set contains volumetric soil water contents and atmospheric boundary conditions (reference evapotranspiration and precipitation) covering almost 12 consecutive months.The meteorological data of the first 45 days were used for model spin-up to allow distortions caused by the calculated initial pressure head distribution to equilibrate.Data from April 2006 to January 2007 were employed for inverse parameter estimation.In total, a number of 1228 daily average water content values were available for the inverse simulations.At the end of the measurement interval between 7 August 2006 and 31 August 2006 the rain gauge was clogged.
To account in the model for the lack in incoming fluxes the amount of water (25 mm) which was removed from the clogged funnel of the device was distributed uniformly over the measured preciptitation values of the time period between 23 August 2006 and 30 August 2006 where the rainguge recorded only 0.4 mm of precipitation while the TDR probe at 0.13 m depth showed a significant increase in soil water content.

Measured data
The bottom section of Fig. 2 displays the volumetric water contents measured at the various depths within the soil profile.As the data set covers a time period of about one year we presume the measured values to be representative for the natural range of water contents for which apparent hydraulic properties need to be estimated.The uppermost TDR probe at 0.13 m depth shows significant fluctuations in volumetric water content during the summer with several wetting and drying cycles.The deeper TDR probes at 0.63 m and 0.92 m depth are not much affected by single rain events.They show a comparably moderate and continuous decrease during the summer with a minimum water content in August www.hydrol-earth-syst-sci.net/13/1953/2009/ Hydrol.Earth Syst.Sci., 13,[1953][1954][1955][1956][1957][1958][1959][1960][1961][1962][1963][1964][1965][1966]2009 and re-wetting in autumn.The water content at 1.16 m depth remains rather constant throughout the year.

Inverse simulations
Inverse simulations were conducted for different model setups.In a first run, the parameters θ r , α, n and K s were estimated for the first four model layers using the preciptitation and reference evapotranspiration values shown in Fig. 2 as upper boundary conditions.This was initially done without invoking a crop coefficient since the grass at our site differs not too much from the reference crop.On average it is slightly shorter and, since we do our measurements for natural boundary conditions, not always well watered.Evapotranspiration was prescribed as given in Sect.2.2.3.As our information on rooting depth (0.08 m) is based on a simple visual inspection from a profile excavation, the influence of rooting depth on the modelled water content evolution was investigated by repeating the simulations using a rooting depth of 0.12 m.The parameter estimates derived from inverse simulations using local gradient-based methods like the trust region algorithm we used in our model study are known to often depend on the choice of the initial parameter set (Vrugt et al., 2008).Hence, to assess the potential dependence of the inversion result on the initial estimate of hydraulic parameters, the simulations for both rooting depths were conducted using sandy loam, loam, silt loam and silt (Carsel and Parrish, 1988) as initial parameter guesses.In total, this leads to eight different inverse simulations.
Independent of the initial parameter guess, all inverse simulations led to similar water content developments.In Fig. 5 measured and simulated values are shown examplary for the sandy loam case.Focusing here on the uppermost TDR probe, measured water contents for both rooting depths clearly cannot be reproduced adequately by the inverse model.While measured and simulated data fit well in autumn and winter and after most wetting events in summer, during April, when evapotranspiration increases (Fig. 2), and again during another drying period in August/September, simulated water contents already show a significant decrease while the measured data still remain at the higher values.Moreover, the model is not able to reach the very low water contents measured during dry periods in summer.As the simulations for the other initial parameter guesses resulted in similar water content phenomenologies, we do not expect the deviation between measured and simulated data to result from the parameter set taken to initialise the model.
As indicated in Sect.2.2.3, one important value affecting the soil water balance in our simulations is the reference evapotranspiration as part of the upper boundary condition.On average, the grass at the test site is shorter than the 0.12 m underlying the FAO Penman-Monteith equation and soil water conditions at our site are not always optimal as it is assumed when calculating j ET 0 .Hence, evapotranspiration may be expected to differ from the calculated values for the reference crop.Furthermore, we trust the rainfall measurements from the weather station and surface run-on on the essentially level experimental field can be neglected.
In a next set of inverse simulations we invoke κ as an additional parameter in the objective function in order to account for differences in j ET 0 compared to the reference crop.In our simple approach, we presume that the general shape of the calculated reference evapotranspiration curve is correct, and that j ET 0 can be scaled with the single crop coefficient κ.Even if the grass at our site is completely developed and generally kept at a certain height, this is a highly simplified model since κ may be expected to vary in time.In general, this may be due to changes in the height of the canopy, reduced evaporation caused by a very dry soil surface etc.By implementing the crop factor in the inversion procedure, our inverse model executes the adjustment of the crop factor automatically with respect to the site-specific conditions.All other parameters were treated as in the previous runs leading to the 17-dimensional parameter vector b={p 1 , ..., p 4 , κ}.Simulations were conducted using four different rooting depths: 0.08 m, 0.10 m, 0.15 m, and 0.20 m in combination with literature values for sandy loam, loam, silt loam and silt (Carsel and Parrish, 1988) for model initialisation.This led to 16 different model realisations (Fig. 6).The goodness of fit of every simulation was evaluated by calculating the root mean square error (RMSE) of all measured and simulated water contents which ranged between 0.013 and 0.038.
Since the different simulations do not result in a common water content phenomenology it is obvious that our model is not able to find a global minimum in the parameter space.We attribute this to the high-dimensionality of the problem and the occurence of local minima as has already been observed in similar modelling studies, e.g., by Ritter et al. (2003) and Hydrol.Earth Syst.Sci., 13,[1953][1954][1955][1956][1957][1958][1959][1960][1961][1962][1963][1964][1965][1966]2009 www.hydrol-earth-syst-sci.net/13/1953/2009/ Fig. 6. Results of the inverse simulations using different initial parameter sets (sandy loam, silt, silt loam, loam (Carsel and Parrish, 1988)) and various rooting depths (0.08 m, 0.10 m, 0.15 m, 0.20 m); black dots: measured water contents (daily average values, only shown for time interval used for inversion); red lines: simulated water contents, RMSE=0.013(n=7); grey lines: simulated water contents, RMSE≥0.014(n=9).Wöhling et al. (2008).However, seven out of 16 runs led to an RMSE of 0.013 and nearly equal temporal variations of water contents at the different TDR probes.Focusing in the following only on these simulations, the inclusion of κ in the objective function improves the results of the inverse simulations significantly and now leads to reasonable fits between measured and simulated data.Particularly, the water contents at the position of the uppermost TDR probe are reproduced much better than in the simulations without κ.
Taking a closer look at Fig. 6 and the corresponding residuals shown in Fig. 7, one still observes significant differences between measured and simulated data.In the uppermost layer, the inverse model is capable of reproducing the measured water contents during wet conditions and also during transitions between wet and dry situations reasonably well.However, after longer dry periods, most prominently at the end of May and beginning of August and a following wetting phase (cf.Fig. 2), the measured water contents are not reproduced well.The observed deviations in our very simple model may originate from a number of factors and we can only speculate on the true reasons.We presume that this discrepancy results from our simple representation of the crop factor which, in fact, should be a function of time.While a reduction of plant transpiration during dry periods is accomplished by the root water uptake function, values for soil evaporation are not adjusted explicitly for all soil water conditions occuring in the profile.During very dry periods, with a dry ground surface, soil evaporation can be expected to tend towards zero.Moreover, evaporation may rise sharply and considerably directly after rain events (Allen et al., 1998).This will potentially prevent an extensive amount of precipitation from infiltrating into the soil profile, predominantly after long, dry periods when the grass leaves are withered and a considerable fraction of the soil surface is exposed to direct solar radiation.For the deeper layers of our model, simulated changes in volumetric water contents are generally  more damped than in the measurements (Fig. 6).A significant difference between measured and simulated water contents occurs in August 2006 at the TDR probes at 0.63 m and 0.92 m depth.During this time the grass at the site had grown significantly higher.This should have resulted in a higher transpiration rate which is not accounted for in our model.At the same depths, the calculated residuals (Fig. 7) show peaks after two distinct wetting events at the beginning and end of October, respectively, when the soil is already pre-wetted after a longer, rainy period during August and September.Here, wetting occus later in the simulations than it does in the measurements which may be an indication for a stronger coupling of the deeper soil sections to the processes in the uppermost layers or, potentially, be a sign for preferential flow.Evapotranspiration is already quite low during this time (cf.Fig. 2).
With respect to the hydraulic parameter estimates inferred from the inverse model, sets with low RMSE values of 0.013 were derived from simulations using sandy loam and loam as initial parameter guesses whereas silt and silt loam parameters lead to higher deviations between measured and simulated data.The estimated hydraulic parameters and values of κ of the runs resulting in an RMSE of 0.013 are listed in Table 3. Plots of the soil water characteristics and hydraulic conductivity functions of the different model layers for these simulations are shown in Fig. 8.All these inverse simulations give reasonable parameter sets for the soil under investigation.However, the inferred parameter estimates are not identical if different initial parameter guesses are used.This is visualised in Fig. 8.The soil water characteristics of model layers 1 and 2 show clear separations within the measured water content range between functions which were derived from simulations initialised with sandy loam or loam as initial parameter guesses.Furthermore, large differences can be observed for the hydraulic conductivity functions of model layers 1 and 4 where the curves vary by about one order of magnitude within the measured water content range.In contrast to the Mualem-van Genuchten parameters, in most cases, rooting depth has a much lower impact on the simulated water content evolutions and soil hydraulic functions of the fine grained soil at our site.This is also manifest in almost identical Mualem-van Genuchten parameters (Table 3) calculated for different rooting depths if the same initial parameter guess is used.This is in accordance to modelling results by Hupet et al. (2002) who found that soil water content has a low sensitivity concerning root water uptake parameters compared to soil hydraulic properties.
The values of κ resulting from the inverse simulations leading to an RMSE of 0.013 are quite similar in all runs and range between 0.59 and 0.62, indicating that the reference evapotranspiration from the FAO Penman-Monteith formula is reduced by the model by about 40%.The reduction of the crop factor appears reasonable since the crop at the Grenzhof Test Site is shorter on average than the reference grass.Midseason values for κ reported by Allen et al. (1998, p. 112) are 0.75 and 0.85 for grazing pasture (extensive grazing) and turf grass (warm season), respectively.These values are valid only for optimal soil water conditions, however.This is not always the case at our site.
The low estimates of the crop factor calculated for our site show the strong influence of κ on the soil water balance.For the time period of our simulations, about 3/4 of the water entering the root zone of the model is immediately removed by transpiration within the uppermost few centimeters of the profile (Fig. 9) and only a small fraction of water is transferred to deeper sections of the soil profile.During summer, water is even taken from greater depths to be available for evapotranspiration.This observation stresses the need for a correct representation of the fluxes passing the soil-plantatmosphere continuum in a model study dealing with natural boundary conditions.
Even though the inverse model is not robust with regard to the resulting Mualem-van Genuchten parameters, the calculated cumulative fluxes across the upper boundary and through the root zone of the model (Fig. 9) from simulations with RMSE=0.013 are very robust throughout all these simulations.The amount of water passing the root zone of the soil profile calculated from these runs ranges between 119 and 138 mm for the simulated time period.Apparently, there exist several different parameter sets that describe the water contents and fluxes equally well.This is a further illustration of Beven's equifinality thesis (e.g.Beven, 1993Beven, , 2006)).

Discussion of overall results
As discussed in the previous section, we explain the observed deviations between measured and modelled water contents by the rather complicated reality which is not represented in detail in our simple model.
Our inverse simulations are based on a traditional local gradient-based method, which does not yield a single global minimum in parameter space.This issue is typically tackled employing global optimisation algorithms like those presented by Duan et al. (1992) or more recently by Vrugt et al. (2008).However, we comment that for the type of problem addressed here -highly simplified representation of complicated reality -statistical and optimisation issues are of lesser importance than analysis of characteristic deviations.
The vegetation model we applied in our study is still fairly simple and neglects the influence of changes in crop height or crop health.Consequently, we still have to expect a considerable uncertainty in the evapotranspiration fluxes calculated by the inverse model.Furthermore, in our case, the conditions are rather simple since the grass is completely developed and generally kept at a certain height, hence we typically do not have to account for different growing stages.However, even in such a simple situation one has to consider the actual evapotranspiration to differ from the one prescibed by the type curve, e.g., due to differences in crop height, albedo, surface resistance or evaporation from the soil surface (Allen et al., 1998).Generally, inverting this simple crop factor is also feasible for other crops or -as an average value -for mixed vegetations.In principle, a temporally changing crop factor could also be obtained from inverting appropriate data that would have to refer to the growing stage of the plants.Unfortunately, such data are not available at our site.Hence, we cannot focus on a more detailed representation of the soil-atmosphere coupling.We still emphasise, however, that when using the FAO Penman-Monteith reference evapotranspiration, the crop factor, also when considered as a constant, should be estimated from the hydraulic state variables.The reason for this is that the average water flux across the upper boundary is of fundamental importance for the soil water balance and the hydraulic state.Any offset, whatever its origin, may lead to characteristic deviations between measured and inverted quantities as is illustrated in Fig. 5.
A further issue that has to be addessed in this context is that soil hydraulic parameters are usually presumed to be static, but may not be so.Before the beginning of the scientific experiments, the Grenzhof Test Site was used for agricultural purposes and the soil was ploughed regularly.Hence, we currently can still expect ongoing textural transitions at least in the uppermost few decimeters of the soil profile.Furthermore, we observe cracks evolving at the soil surface during dry conditions which also influence the hydraulic properties of the soil at our site.Although these processes are not included in our model study, we have to consider their influence on the temporal development of observed soil water contents.

Extension to the field scale
In this study, estimation of the hydraulic properties of different soil layers was restricted to one single soil profile, which does not automatically lead to a parameterisation which is representative for an entire field with extents of hundreds of meteres to kilometres.The well-known heterogeneity of soils impedes this.However, the estimated values are relevant at a scale that is appropriate to simulate water flow at the field scale.The correct approach now would be to run many such profiles in order to obtain the parameter field with the desired extent and resolution, and to run a corresponding high-resolution, three-dimensional numerical simulation.In our understanding, the sometimes proposed "field-averaged soil hydraulic material properties" is -for a heterogeneous field -an invalid concept since the very definition of material properties demands local equilibrium at the scale of representation.This, in general, cannot be ascertained for the Richards equation under typical atmospheric forcing (Roth, 2008).

Summary and conclusions
We estimated the hydraulic properties of a layered field soil from numerical inversion using a time series of in situ measured volumetric soil water contents.The model was driven by the transient fluxes across the soil-atmosphere-interface which were measured by an automatic weather station located next to the instrumented soil profile.The inverse estimation of hydraulic parameters in combination with a crop factor which scales the FAO Penman-Monteith reference evapotranspiration, while based on a rather simple process model, results in quite a reasonable agreement between measured and simulated data.
One important observation revealed by our simulations is the importance of the fluxes passing the soil-vegetationatmosphere continuum.Here, rather accurate values of the time-dependent evapotranspiration fluxes are necessary to better represent the upper boundary of the model which drives the simulation.This is so far not considered in many inverse modelling studies and we suggest to include the adjustment of the evapotranspiration fluxes in the inverse parameter estimation.
Besides of being very cost effective, the method circumvents two crucial issues: (i) extraction of undisturbed soil samples and (ii) estimating field-scale properties from those determined in the laboratory.Similarly comprehensive data sets are recorded in many transient field experiments by many monitoring stations, hence the method is broadly applicable.If this is not the case, installation of sensors is quick, and cheap if manual reading is chosen which suffices for many situations.The method has clear limitations as was illustrated by the postulated preferential flow.This is more in the model formulation, however, and not so much in the actual method.While further studies will have to substantiate the claim, we see the prospect of abandoning laboratory experiments entirely as far as field-scale understanding and prediction is the goal.
Based on this study, we propose that for a site where continuous measurements of soil water content and meteorological data are available and the soil is not subject to complex processes like, e.g., preferential flow, inverse estimation of apparent hydraulic parameters from in situ measurements is a good alternative to traditional laboratory methods.

Fig. 2 .
Fig. 2. Daily average values of reference evapotranspiration, daily sums of precipitation and daily average values of volumetric soil water content measured at the Grenzhof Test Site between Febuary 2006 and January 2007.Meteorological data from the first 1.5 months (Febuary to March 2006) were used for model spin-up, data starting from April 2006 was taken for inverse parameter estimation.Values for precipitation between 23 August 2006 and 30 August 2006 were corrected in order to account for the lack in incoming fluxes caused by the clogged rain gauge (Sect.2.2.4).

Fig. 5 .
Fig. 5. Results of inverse simulations for the temporal water content evolution at 0.13 m depth using measured precipitation and calculated reference evapotranspiration (j ET 0 ) data to represent upper boundary fluxes.Results calculated for two different rooting depths are compared to the measured water contents.During summer (high j ET 0 ) the model is not able to reproduce the evolution of the measured data (only shown for time interval used for inverse parameter estimation).

Fig. 7 .
Fig. 7. Residual water contents (measurement-simulation) at the different TDR probes calculated for runs with RMSE between measured and simulated water contents of 0.013.For comparison, measured water contents are displayed in the topmost frame of the figure.

Fig. 8 .
Fig. 8. Soil-water characteristics and hydraulic conductivity functions derived from inverse simulations with RMSE=0.013.Grey bars indicate the range of measured water contents, various lines for both soil textures indicate simulations with different rooting depths.

Fig. 9 .
Fig. 9. Cumulative fluxes calculated for the time period of inversion, top: surface flux (precipitation minus evaporation from the soil surface), middle: root water uptake, bottom: water flux passing the root zone; red: RMSE=0.013;grey: RMSE≥0.014.Considering only the simulations with RMSE of 0.013, during the time period under investigation approximately 3/4 of the water infiltrating the soil profile is removed by root water uptake.

Table 1 .
Configuration of weather station: sensors used for this study.Rotronic AG, Switzerland; b R. M. Young, Traverse City, MI; c Kipp & Zonen, Delft, The Netherlands; d Campbell Scientific, Logan,

UT Fig. 1. Soil profile from the Grenzhof Test Site (left). The different colours on the scale indicate intervals of
• C] and u 2 is the wind speed [m s −1 ], both at 2 m height, e s is the saturation vapour pressure [kPa], e a is the actual vapour pressure [kPa], e s −e a is the saturation vapour pressure deficit [kPa],

Table 3 .
Apparent hydraulic parameters and crop coefficients inferred from simulations with RMSE=0.013.Values in brackets indicate 95% confidence intervals.