Repositorio Institucional de Documentos

Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es  Trabajo Fin de Máster Estimation of water content profiles by inverse analysis of TDR waveforms: application to infer soil hydraulic properties Autor/es Carolina Peña Sancho Director/es Dr. David Moret-Fernández Dr. César González-Cebollada Escuela de Ingeniería y Arquitectura 2013 Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 1 En primer lugar, quisiera agradecer al Dr. Borja Latorre, de la Estación Experimental Aula Dei, por su inestimable ayuda en todo lo referente a programación de este trabajo. También, expresar mi agradecimiento al Dr. Francisco Lera, de la Universidad de Zaragoza, porque sin su colaboración este trabajo nunca hubiera salido adelante. Por último, agradecer también la ayuda y compañía de Mariví, de Pepa, y de Ana, de la Estación Experimental Aula Dei. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 2 Abstract The measure of soil hydraulic properties, water retention curve (θ(ψ)) and hydraulic conductivity (K) results of paramount importance for hydrological processes simulation. The Time Domain Reflectometry (TDR) technique is a worldwide used technique that allows non destructive measurements of both soil volumetric water content (θ) and electric conductivity (σ). The goal of this work is to develop a new TDR based methodology to estimate soil hydraulic parameters (α, n and K) by inverse analysis of WCP’s dynamics under water falling infiltration experiments. The WCPs are estimated from the inverse analysis of the TDR waveforms using a physical electromagnetic propagation model. The α, n and K are subsequently calculated using a HYDRUS-1D-Matlab interface by inverse analysis of estimated TDR-WCPs. This interface calculates the hydraulic parameters for the best fitting between the recorded TDR-WCP set and those simulated by HYDRUS-1D. To this end, a brute-force optimization method is employed, which allows sweeping a wide range of hydraulic parameters. The method was tested on three different porous media (2-mm sieved loam soil, sand, and glass microspheres) during a water falling infiltration process. The hydraulic properties estimated with this method were compared to those measured for the same porous media using conventional laboratory methods: TDR-pressure cell and mini-disc infiltrometer. Although satisfactory estimations of K were obtained, inaccurate n and α value were observed. These discrepancies could be attributed (i) the unimodal instead of a bimodal function used by HYDRUS-1D; (ii) the soil hysteresis phenomena; (iii) uncertainties on the “bridge” function used to estimate the WCP from the modeled TDR waveforms. New efforts are needed from improve the accuracy of this method and testing it in undisturbed soil cores. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 3 Resumen La medida de las propiedades hidráulicas del suelo, curva de retención (θ(ψ)) y conductividad hidráulica (K) tiene una importancia fundamental para la simulación de procesos hidrológicos. La técnica de Reflectometría de Dominio Temporal (TDR) es una herramienta ampliamente utilizada para la medida no destructiva de contenido volumétrico de agua en el suelo (θ) y de conductividad eléctrica (σ). El objetivo de este trabajo es desarrollar una nueva metodología basada en el uso de la técnica TDR para estimar las propiedades hidráulicas del suelo (α, n y K) por análisis inverso de la dinámica de los perfiles de humedad (WCPs ) durante un proceso de infiltración de agua. Los WCPs se estiman a partir del análisis inverso de ondas TDR empleando un modelo físico de propagación electromagnética. Posteriormente, los parámetros α, n y K se calculan empleando una interfaz HYDRUS-1D-Matlab por medio del análisis inverso de los TDR-WCPs. Esta interfaz calcula los parámetros hidráulicos a partir del mejor ajuste entre los TDR-WCP registrados y los simulados por HYDRUS-1D. Para este fin, se emplea el método de optimización de fuerza bruta, el cual permite barrer un rango amplio de parámetros hidráulicos. El método fue probado en tres medios porosos distintos (tierra franca tamizada a 2 mm, arena y microesferas de vidrio) durante un proceso de infiltración. Las propiedades hidráulicas estimadas con este método se compararon con aquellas medidas en los mismos medios porosos empleando técnicas convencionales de laboratorio: cámaras de presión-TDR y mini-infiltrómetro de disco. Aunque se obtuvieron resultados satisfactorios para la medida de K, los resultados obtenidos para la medida de n y α fueron imprecisos. Estas discrepancias se pueden atribuir a las siguientes causas (i) el uso de una función unimodal en lugar de una bimodal en HYDRUS-1D; (ii) el fenómeno de histéresis del suelo; (iii) incertidumbres en la función “puente” empleada para estimar WCP a partir del modelo de simulación de ondas TDR. Este método necesita nuevos esfuerzos para mejorar su precisión y así poder probarse en muestras de suelo inalterado. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 4 TableofContents 1. Introduction..................................................................................................................................6 2. Theory..........................................................................................................................................9 1. Soil water flow....................................................................................................................9 1.1. Soil hydraulic properties ...................................................................................................9 1.2. Soil water infiltration ......................................................................................................10 1.3. Subsurface soil water flow..............................................................................................12 1.4. Hydrus 1D Software Package version 4.15.....................................................................14 2. Time Domain Reflectometry (TDR) waveforms analysis to estimate the soil water content and the bulk electrical conductivity...............................................................................16 2.1 Estimations of volumetric water content and bulk electrical conductivity with the graphical method....................................................................................................................16 2.2 Numerical model to estimate the volumetric water content.............................................17 3. Optimization techniques and sensitivity analysis..............................................................20 3.1. Optimization techniques..................................................................................................20 3.2. Sensitivity analysis..........................................................................................................24 2. Material and Methods ................................................................................................................25 2.1 Development of a TDR based method for water content profiles estimation......................25 2.2 Matlab interface to HYDRUS 1D........................................................................................27 2.2.1 Hydrus 1D - Matlab interface.......................................................................................27 2.2.2 Optimization process for hydraulic parameter estimation.............................................31 2.2.3 Zoom optimization........................................................................................................32 2.2.4 Sensitivity analysis........................................................................................................33 2.3 Validation of the new method for estimate soil hydraulic properties ..................................33 2.3.1 Measurement of soil hydraulic properties.....................................................................34 2.3.2 Experimental design to estimate soil hydraulic properties by inverse analysis of the water content profiles.............................................................................................................37 3. Results and discussion................................................................................................................38 Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 5 3.1 Experimental results.............................................................................................................38 3.1.1 Particle size distribution of the porous media ...............................................................38 3.1.2 Measurement and modeled water retention curve.......................................................39 3.1.3 Cumulative infiltration curve and soil hydraulic conductivity.............................42 3.2 Water content profiles dynamics modeling..........................................................................44 3.2.1 Measured and modeled TDR waveforms during falling infiltration experiments.........44 3.1.2 Dynamics of the water content profiles during the falling infiltration experiment measured by TDR and modeled by Hydrus 1D .....................................................................46 3.2 Matlab-HYDRUS 1D optimization interface ...............................................................47 3.2.1 Estimation of hydraulic parameters (α, n, K) for porous media....................................47 3.2.2 Sensitivity analysis of hydraulic parameters........................................................50 4. Conclusions................................................................................................................................54 References......................................................................................................................................55 Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 6 1.Introduction  The determination of the soil hydraulic properties is of paramount importance in many scientific fields such as agronomy, hydrology and environmental science. The water flow into the soil depends on its ability to transmit water through the porous medium. This is function on the pore-size distribution, tortuosity, shape and degree of interconnection of the water-conducting pores in the porous medium. The parameters that define the water flow into the soil are the hydraulic conductivity, K and the water retention curve  (ψ) (Dane and Hopmans, 2002). The water retention curve θ(ψ) is the relationship between the soil volumetric water content (θ) [m3/m3] and the matric potential (ψ) [KPa]. The shape of this function depends on the soil aggregates and particle size distribution. As suggested by Guérif et al. (2001) the soil porosity can be considered as (i) textural porosity that occurs between the primary mineral particles and depends on organic matter content and soil texture (Dane and Hopmans, 2002), and (ii) structural porosity, sensitive to soil management factors and comprised by microcracks, cracks, bio-pores, and macrostructures produced by tillage (Dexter, 2004). Estimates of θ(ψ) require pairs of ψ and θ measurements. The most common laboratory technique for estimating θ(ψ) is the pressure plate extractor. In this method, the water retention information is obtained by bringing the soil sample to equilibrium by applying a constant pressure gradient across the soil, driving water movement while preventing air entering into a pressurized chamber (Dane and Hopmans, 2002). The θ is commonly calculated from the measured soil gravimetric water content and dry bulk density. Although the pressure plate extractor can house undisturbed soil contained in a metallic cores (Wraith and Or, 2001), this technique is mainly applied on sieved soil samples placed on 2 cm high rubber cylinders. Soil saving makes that pore size-distribution of soils and consequently θ(ψ) substantial changes regarding to the undisturbed field conditions (Moret-Fernández er al. 2012). However, water retention curves of undisturbed soil samples are highly desirable to more accurate modeling of soil water flow and water balances. On the other hand, the gravimetric method used to measure θ is trying, tedious and time consuming. To solve this limitations, Moret-Fernández et al. (2012) developed a Time Domain Reflectometry (TDR) based pressure cell to determine θ(ψ). This consists of a 50-mm internal diameter stainless steel cylinder attached to a porous ceramic disc, closed at the ends with two aluminum lids, and longitudinally crossed by a stainless steel rod. Although this method worked with undisturbed soil samples and allowed simplifying the water content measurement, the discontinuous sampling of this method makes the estimate of a representative soil water retention curve to be time consuming (up to 14 days per retention curve). Soil hydraulic conductivity (K) is a measure of the soil ability to transmit water when soil is submitted to a hydraulic gradient. K is a function of soil water content, the hydraulic head, and the Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 7 flux across the upper boundary of a soil compartment (Dane and Hopmans, 2002). Sorptivity (S) is a measure of the ability of an unsaturated porous medium to absorb or store water as a result of capillarity (Philip, 1957; Dane and Hopmans, 2002). So far, different laboratory and field methods to estimate K on both disturbed and undisturbed soil samples have been developed. For instance, some standard laboratory methods are the constant head soil core tank method, the falling head soil core tank method, or the steady flow soil column method (Dane and Hopmans, 2002). Field methods, which allow in situ determination of K and S, are mainly based on the infiltrometry technique. Field methods cover from the simplest single or double ring methods to the more complex Guepth permeameter or the tension disc infiltrometers. Over the last two decades tension disc infiltrometers (Perroux and White, 1988) have become very popular devices for in-situ estimates saturated an unsaturated K and S (White et al., 1992) and macropore flow contribution (Angulo-Jaramillo et al., 2000). Tension disc infiltrometers consists of a disc base covered by a membrane, a graduated water-supply reservoir and a bubble tower with a moveable air-entry tube that imposes the pressure head at the cloth base (Perroux and White, 1988). The hydraulic properties (K and S) are commonly calculated from the measured cumulative infiltration curve. Two different methods are so far available: the steady-state and the transient water flow methods. Compared to the standard the steady-state water flow method (Ankeny et al., 1991), the transient water flow procedure, that requires shorter experiments, involves smaller sampled soil volumes and consequently more homogeneous and initial water uniformity (Angulo-Jaramillo et al., 2000). Several simple expressions have been developed to estimate the soil hydraulic parameters from the transient water flow (Warrick and Lomen, 1976; Warrick, 1992, Vandervaere et al., 2000; Zhang, 1998). However, based on the quasi-exact analytical form of the 3D cumulative infiltration curve from the disc infiltrometer (Haverkamp et al. 1994), Latorre et al. (2013) proposed a new method to calculate the K and S from the numerical solution of the complete Haverkamp et al. (1994) model. This new procedure, which results robust enough, allowed better estimates of the hydraulic properties. Simulation models are interesting to simulate the water balance in soil-crop systems (Connolly, 1998). Over the large amount of soil physical models so far available, the HYDRUS is a worldwide used model that numerical solves the Richards flow equation in a porous media (Simunek, 2009). The HYDRUS-1D code may be used to analyze water and solute movement in unsaturated, partially saturated, or fully saturated porous media. The flow region itself may be composed of nonuniform soils. Flow and transport can occur in the vertical, horizontal, or in a generally inclined direction. The water flow part of the model considers prescribed head and flux boundaries, as well as boundaries controlled by atmospheric conditions, free drainage, or flow to horizontal drains (Simunek, 2009). Time Domain Reflectometry (TDR) has become a worldwide standard technique that allows simultaneous, accurate, and non-destructive estimations of volumetric water content and bulk Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 8 electrical conductivity (  a) (Topp and Ferré, 2002). The TDR instrument launches an electromagnetic pulse along a probe embedded in a porous medium, and the signal is displayed as a TDR waveform in which the voltage (V) or reflection coefficient (  ) is expressed as a function of time (t). While the travel time along the probe depends on the probe length and the apparent permittivity in the vicinity of the probe, the value of which is highly correlated to  , the  is also dependent on the  surrounding the TDR probe (Topp and Ferré, 2002). Two different approaches to determine both  and  from the TDR waveform are currently available. The first, which estimates the soil parameters by a graphical analysis of the TDR waveform, uses the travel time (tL) to estimate the dielectric permittivity, and the attenuation of  to assess the  a. The second and more sophisticated method to estimate  and  is based on the modeling of the TDR waveform using the physical properties of the system (Friel and Or, 1999; Heimovaara et al., 2004; Greco, 2006). In this approach,  and  are modeled parameters estimated by fitting them against measured waveforms. Although the method requires significant computer resources, the number of outliers due to erroneous analyses is substantially reduced (Heimovaara et al., 2004). Optimization is central to any problem involving decision making, whether in engineering or in economics. The task of decision making entails choosing among various alternatives. This choice is governed by the desire to make the ‘best’ decision. The measure of goodness of the alternatives is described by an objective function or performance index (Chong E.K.P & Zak S.H. 2013). In the simplest case, an optimization problem consists of maximizing or minimizing a real function by systematically choosing the input values from within an allowed set and computing the value of the function. The generalization of optimization theory and techniques to other formulations comprises a large area of applied mathematics. More generally, optimization includes finding “best available” values of some objective function given a defined domain, including variety of different types of objective function and different types of domains (Chong E.K.P & Zak S.H. 2013). The objective of this work is to develop a new method to estimate the hydraulic parameters of soil from the inverse analysis of the water content profiles (WCPs) under water falling infiltration experiments. The WCPs were obtained from the inverse analysis of the TDR waveforms using a physical electromagnetic propagation model. An interface HYDRUS-1D-Matlab was developed to estimate hydraulic parameters from TDR measured water content, and the method was tested on three different porous media.   Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 15 gradients in the water retention function are too big, and numerical convergence problems in HYDRUS 1D simulations. For the same reason, the initial water content in the “saturated zone” e.g. on the top of the column for an infiltration process, must be smaller than the water content at saturation own of the porous media. The governing flow and transport equations are solved numerically using standard Galerkintype linear finite element schemes, or modification thereof. The program is a one-dimensional version of the HYDRUS-2D and HYDRUS (2D/3D) codes simulating water, heat and solute movement in two or three-dimensional variably saturated media (Šimůnek et al., 1999; 2006a,b), while incorporating various features of earlier related codes such as SUMATRA (van Genuchten, 1978), WORM (van Genuchten, 1987), HYDRUS 3.0 (Kool and van Genuchten, 1991), SWMI (Vogel, 1990), SWMI_ST (Šimůnek, 1993), HYDRUS 5.0 (Vogel et al., 1996), and HYDRUS- 1D, version 3.0 (Šimůnek et al., 2005,Simunek, 2012). In addition, HYDRUS 1D implements a Marquardt-Levenberg type parameter estimation technique for inverse estimation of soil hydraulic and/or solute transport and reaction parameters from measured transient or steady-state floe and/or transport data (Simunek, 2012). Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 16 2. TimeDomainReflectometry(TDR)waveformsanalysistoestimate thesoilwatercontentandthebulkelectricalconductivity  2.1 Estimations of volumetric water content and bulk electrical conductivity with the graphical method The transit time of the TDR pulse propagating one return trip in a transmission line of length L (m), tL, is expressed by: c L tL a ε2  (22) where c is the velocity of light in free space (3 x 108 m s-1) and  a is the apparent permittivity of the medium (Topp and Ferré, 2002). The tL value is calculated as the distance between the time at which the signal enters the TDR rods (first peak) and the time when the trace arrives at the end of the TDR probe, also denoted second reflection point or end point. These points can be manually determined or calculated using a computer algorithm to find the end point. In this case, the most used procedure is the “tangent method” (Heimovaara, 1993). The volumetric water content  , can be calculated from  a according to Malicki equation (1996).    18.117.7 159.0168.0819.0 ),( 2 **   a (23) where ρ is the soil bulk density (Malicki, 1996). The soil bulk electrical conductivity (  a) estimated with the graphical long-time TDR waveform analysis is calculated according to (Giese and Tiemann, 1975):              Scale, Scale, aρ1 ρ1 σ r p Z K (24) where Zr is the output impedance of the TDR cable tester (50 Ω), Kp (m-1) is the probe-geometry- dependent cell constant value, and Scale, ρ is the scaled steady-state reflection coefficient corresponding to the ideal condition in which there is no instrument error or cable resistance. The Scale, ρ is calculated using the equation described by Lin et al. (2008): 1 )1)(())(1( ))(( 2 ,      airSCairairSC airSCair Scale       (25) Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 17 where  , air  and SC  are the long-time reflection coefficient measured in the studied medium, in the air and in a short-circuited probe, respectively. The reflection coefficient  , as a function of time, t, is defined as    i VV VtV t   0 0 ρ-1    +1 (26) where V(t) is the measured voltage at time t, V0 is the voltage in the cable just prior to the insertion of the probe (standard impedance value of 50 ), and Vi is the incident voltage of the cable tester prior to the pulse rise. 2.2 Numerical model to estimate the volumetric water content The TDR signal ρ(t) is the transient response of the cable-probe-soil set to the cable tester excitation signal. The cable and probe will be modelled as lossy transmission lines in the frequency domain. Fourier analysis (FFT) will be used (Heimovaara, 1994; Heimovaara, 2004; Huebner and Kupfer, 2007) with direct and inverse FFT algorithms for switching from time to frequency domain and vice-versa. The excitation signal used in the modelling process is the actual cable tester output measured in open circuit. The frequency domain transfer function of the soil-probe-cable set is that of a voltage divider constituted by the output impedance of the cable tester Zr (nominally 50 Ω) and the frequency-dependent input impedance of the cable-probe-soil set, Zi. Four distributed parameters are used to characterize transmission lines (Ramo et al., 1984): capacitance C (F m-1), inductance L (H m-1), conductance G (S m-1) and resistance, R (Ωm-1). Due to geometrical considerations, for lines of uniform cross section and for linear media:   LC    //  CG (27) The characteristic impedance, Zo (Ω), and the propagation constant ϒ(m-1) at angular frequency ω are then obtained as:    CjGLjRj CjG LjR Z        0(28) where α is the attenuation constant (Np m-1) and β is the phase constant (rad m-1). For ideal lossless lines, R=0 and G=0 and thus:    jLCjj C L Z 0(29) Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 18 where ν is the phase velocity. The input impedance of a transmission line with known characteristics (Zo and ϒ) of length l connected to a load impedance ZL is calculated as: l l L L i   tanhZZ tanhZZ ZZ 0 0 0  (30) The input impedance of the cable-probe-soil set Zi, is computed in a two-step process. Firstly, we apply the equation to obtain Zp, the input impedance of the probe inserted in the soil as a transmission line of length l=lp ending in an open circuit load (ZL->infinity), using the probe’s characteristic Zop, and  p values corresponding to a given  and  apair. Secondly, Zi is obtained, again using the equation to compute the input impedance of the coaxial cable as a transmission line of length l=lc ending now in a load impedance ZL= Zp, using the coaxial cable’s characteristic Zoc, and  c  values. We have used a coaxial cable of type RG58, with nominal Zo = 50  and   = 0.66 c. Using the equation we obtain Lc = 250 nH m-1 and Cc = 100 pF m-1. In the TDR frequency range, skin effect losses are the dominant ones and give rise to a series resistance, Rc, and to an extra external inductance, Lc2. Both terms are frequency-dependent (Nahman 1972). The values obtained from best fits to TDR measurements of the coaxial cable ending in open and short circuit are:  mHLmmR cc / 177 /177.040 2     (31) The transmission line parameters for lossless three-rod probes in air–three identical cylindrical rods of length lp, radius b and center-to-center spacing s – have been derived from the calculations of Ball (2002). The characteristic impedance in a vacuum or air, Zp0, is very well approximated by the following expression, where d = b/s.        3 0 0 02 1 ln 4 1 d Zp    (32) Cp0 and Lp0, can be derived from their respective equations, for the probe in air as     3 0 0 3 0 02/1ln 42/1ln 4dL d Cpp     (33)  Skin effect losses can be neglected for our short probe lengths. Short probes need a correction of their actual length to an effective, longer one, due to the fringing of the electromagnetic field at the probe’s open end. We include this correction, adding an extra length, double that estimated by Green and Cashman (1986) for two-rod probes. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 19 When the probe is inserted in a lossy soil with bulk conductivity  , G is obtained substituting 0 with  . A direct estimation of  by the Castiglione and Shouse (2003) method requires the cell constant value Kp, usually obtained by a calibration procedure. As G is now known, the theoretical expression for the cell constant of trifilar probes can be used instead. Here LM is the effective length of the probe, longer than the actual length:        3 2 1 4 1 d ln LπGl σ K Meff p(34) Dielectric effects, including losses, are incorporated by substituting  0 in the above equation, with the soil complex permittivity  c=  ’–j  ’’To estimate  cwe first compute the frequencydependent complex permittivity of pure water  w(  ) at a given temperature, following Meissner and Wentz (2004). For a given  we obtain  a (  ) with the Malicki equation and finally:  0 01 0 0aw aa aa ac          (35) where  a0=  a (  =0) and  a1=  a (  =1).    Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 20 3. Optimizationtechniquesandsensitivityanalysis 3.1. Optimization techniques Optimization is central to any problem involving decision making, whether in engineering or in economics. The task of decision making entails choosing among various alternatives. This choice is governed by the desire to make the ‘best’ decision. The measure of goodness of the alternatives is described by an objective function or performance index (Chong E.K.P & Zak S.H. 2013). In the simplest case, an optimization problem consists of maximizing or minimizing a real function by systematically choosing the input values from within an allowed set and computing the value of the function. The generalization of optimization theory and techniques to other formulations comprises a large area of applied mathematics. More generally, optimization includes finding “best available” values of some objective function given a defined domain, including variety of different types of objective function and different types of domains (Chong E.K.P & Zak S.H. 2013). An optimization problem can be represented in the following way: - Given a function f : A → R from a set A to the real number - Sought: an element X0 in A such that f(X0) ≤ f(X) for all x in A, this is called ‘minimization’; or such that f(X0) ≥ f(X) for all x in A, this is called ‘maximization’. (36) Such a formulation is called an optimization problem, many real-world and theoretical problems may be modelled in this general framework (Chong E.K.P & Zak S.H. 2013). Typically, A is some subset of the Euclidean space Rn, often specified by a set of constraints, equalities or inequalities that members of A have to satisfy. The domain A of f is called the ‘search space’, while the elements of A are called ‘candidate solutions’ or ‘feasible solutions’. The function f is often called an ‘objective function’, and is feasible solution that minimizes, or maximizes the solution (Chong E.K.P & Zak S.H. 2013). By convention, the standard form of an optimization problem is stated in terms of minimization. Generally, unless both the objective function and the feasible region are convex in a minimization problem, there may be several local minima where a local minimum x is defined as a point for which there exists some δ > 0. So that for all x such that ||x –x*|| ≤δ (37) the expression Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 21 f(x*) ≤ f(x) (38) holds, that means, on some region around x* all of the function values are greater than, or equal to the value at that point. Local maxima are defined similarly (Chong E.K.P & Zak S.H. 2013). A large number of algorithms proposed for solving non-convex problems are not capable of making a distinction between local optimal solutions and rigorous optimal solutions. The branch of applied mathematics and numerical analysis that is concerned with the development of deterministic algorithms that are capable of guaranteeing convergence in finite time to the actual optimal solution of a non-convex problem is called, global optimization (Chong E.K.P and Zak S.H. 2013). Due to the optimization problem of this work is a non-convex problem, an unconstrained minimization problem will be solved. The special characteristic of this problem is that the solution of the three elements vector X need not satisfy any constraint. Several methods are available for solving an unconstrained minimization problem. These methods can be classified into two broad categories: as direct search methods and descent methods. The direct search methods require only objective function evaluations and do not use any partial derivatives of the function in finding the minimum and hence are often called nongradient methods. These methods are most suitable for simple problems involving relatively small number of variables, and they are, in general, less efficient than the descent methods. The descent techniques require, in addition to function evaluations, the evaluation of a first and possibly higher order derivatives of the objective function, these techniques are also known as gradient methods. Some of the direct search methods most commonly used are: brute force method, random search method, univariate method, pattern search method and Simplex method. And some of the descent methods most commonly used are: steepest descent method, conjugate gradient method, Newton’s method, and Gauss-Newton method (Rao S.S., 1984). There are also a method called Marquardt-Levenberg method, that is actually a combination of two minimization methods: the gradient descent method and the Gauss-Newton method. The objective function to be minimized in this work is a three parameter vector.                       s K n a x x x X 3 2 1 (39) In this case, the function to be minimized is a nonlinear function of error between WCP set estimated by inverse analysis of the TDR waveforms and WCP dynamics calculated by HYDRUS 1D simulations. The function of error can be calculated by the Root Mean Square Error (RMSE) that is a difference between values predicted by the model or estimator and the values actually observed. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 22 n yy RMSE n t tt    1 ^2 )( (40) These individual differences are called residuals when the calculations are performed over the data sample, and are called ‘prediction errors’ when computed out-of-sample. The RMSE serves to aggregate the magnitudes of the errors in predictions for various times into a single measure of predictive power and is a good measure of accuracy to compare forecasting errors of different models for a particular variable. As mentioned, the adequate imposition of initial conditions in water contents is of major importance for HYDRUS 1D simulations. There is not existing a unique way to calculate the suitable initial conditions. To solve this problem, the Midpoint rule has been employed. This consists of a root-finding method that repeatedly bisects an interval and selects a subinterval in which the root must be lie for further processing. The method calculates a middle point in the considered interval xr as: 2 31 2 xx x   (49) According to Bolzano’s theorem, the method is applicable when the equation f(x)=0 must to be solved, where f is a continuous function defined on the [x1,x3] interval and f(x1) and f(x3) have opposite signs and bracket a root. In this case x1 and x3 are said to bracket a root since, by the intermediate value theorem, the f is a binary function that indicates if HYDRUS 1D works (it gives a positive error), or stops (it gives a negative error). 1 e1 e2 e3 -1 Figure 1. Midpoint rule for the function of error Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 23  In order to minimize the function of error obtained, the ‘brute force’ search method will be e employed. This is an algorithm that tries all possible solutions until it finds an acceptable one or until a pre-set maximum number attempts. A brute-force optimization algorithm evaluates value after value for a given loop, and return the value with the optimal result. In this case, the optimal result is the minimum of the given function of error. End xf xf If )(min min)(   (50) Although the number of possible states of the system increases exponentially with the number of dimensions, brute force methods have the benefits that they are simple to implement, and in the case of discrete systems, all possible states are checked. As consequence, brute-force methods are useful to determine the general shape of the function to minimize, and for realizing the rear sensitivity analysis. Other optimization methods such as multivariate Newton-Rhapson method, or the Marquardt- Levenberg algorithm could be used. Other methods based on the search of the maximum gradient of a given multidimensional function can be used for accelerate the optimization process (Eq. 51). For instance, we have, the steepest descend method, which searches the negative of the gradient vector as a direction for minimization; the conjugate gradient method, which uses the quadratic convergence for search the minimum of the function (Rao, 1984); or the Marquardt-Levenberg algorithm, that is a combination of the gradient descent and the Gauss-Newton method (Gavin 2011). This last method is the optimization technique employed by HYDRUS-1D. 0)(':  ii xfx (51) However, although these methods have the advantage of their fast time calculation, they have the disadvantage that they can may find a local minima, so they are not be able to find the global minimum, this is the reason that led us to use the brute force method for establish the function of error. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 24 3.2. Sensitivity analysis In most of the practical problems, we are interested not only in the optimal solutions of the optimization problem, but also how the solutions changes when the parameters change. The study of the effect of discrete parameters changes on the optimal solution is called the sensitivity analysis (Rao, 1984). One way to solve the effects of changing parameters is solving a series of new problems. In general, when a parameter is changed, it results in one of the three cases (Rao, 1984):  The optimal solution remains unchanged.  The basic variables remain the same but their values are changed.  The basic variables as well as their values are changed. In this work, a sensitivity analysis will be made by performing all the optimization process sweeping the range of the involved parameters one by one. The final goal of the sensitivity analysis is to perform a graph with the effect of each discrete parameter change on the value of the function of error.  Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 31 2.2.2 Optimization process for hydraulic parameter estimation The HYDRUS 1D-Matlab interface also includes the optimization process to estimate hydraulic parameters of porous media. The optimization algorithm searches, using the bruteforce method, the minimum of the function of the calculated RMSE for each α, n, and Ks parameter combination. The goal of this brute-force method is to create a map of error in order to analyze, separately, the shape of the functions of error for each parameter, in order to subsequently consider the possibility of employ some of the Optimization Toolbox of Matlab 2010. The objective of these tools is accelerating the optimization process. Firstly the optimization process imposes a wide range for minimum and maximum values of α, n, and Ks variables. It should be noted that, due to logarithmical nature of α and Ks variables, the code runs the logarithm of these variables during the optimization process. αi αf n i n f K i K f log(0.0001) log(0.0125) 1.0001 log(5.84) log(0.0002) log(0.01) Table 1. Parameter’s range swept during the brute-force optimization process. The number of attempts of the brute-force method is introduced, so the size of the mesh is calculated.         20 20 20 n n n K n  (62)            nif nif nif KKKdK nnndn d / / /  (63) Then, a loop with α, n, and Ks variables is established for calling the HYDRUS 1D simulation. This function calculates also the RMSE between the WCP set done by each HYDRUS 1D simulation of the loop and the WCP set coming from the TDR inverse analysis model. The RMSE minimized with the brute-force method, allows obtaining a complete map of the function of error for the swept parameters range. The number of iterations corresponding to the minimum of α, n, and Ks respectively is also saved. For i1 = 0:αn Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 32    di i* 1  For i2 = 0:nn dninn i* 2   For i3 = 0:Kn dKiKK i* 2   error =simulation1[exp(α) n exp(Ks) θs θr ] If error < min min = error min(α) =α min(n) =n min(Ks) = Ks miniα = iα minin = in minik = ik (64) Finally, the code plots in the same graphic the dynamics of the WCP simulated by HYDRUS 1D and those obtained from the inverse analysis of the TDR waveforms. The graph changes when reducing RMSE for the two models is observed. 2.2.3 Zoom optimization Once the minimum error, the parameter values corresponding to this minimum and the number of iterations for each parameter corresponding to the minimum are obtained, a “zoom optimization” is applied. This allows finding the global minimum, and consequently preventing possible local minima in the vicinity of the global minimum. The zoom optimization code is very similar to the global optimization code. The difference lies in the calculation of the boundary values of α, n, and Ks parameters for the optimization process. The iteration corresponding to the minimum value of each parameter is added to the initial parameter value and then, a few steps are added, subtracted and multiplied by the global size of the mesh. This allows establishing the “optimization window” around the global minimum. For instance, the process to calculate the boundary values and the size of the mesh for α parameter runs as follows αf= αi + (miniα + 4)*dα αi= αi + (miniα - 4)*dα (65) Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 33 It is important to note that the size of the mesh must be recalculated for the zoom optimization with the new boundary values. Then, the optimization process is similar to the global optimization process. 2.2.4 Sensitivity analysis The sensitivity analysis requests that an optimization loop for sweep each hydraulic parameter range to be separately performed to execute the HYDRUS 1D simulations. For each loop, just one parameter is swept on their wide range, and the other two parameters are fixed on their minimum values coming from the zoom optimization process. For example, the loop used for the α parameter sensitivity analysis runs as follows α n = 200 dα = (αf – αi)/ α n αmin = a nmin = b Kmin = c For iα = 0: α n α = αi + iα* dα error =simulation1[exp(α) n exp(Ks) θs θr ] End (66) Firstly the attempts number is established, next the size of the mesh calculated, the minimum values for each parameters fixed, and finally the optimization process for α parameter is implemented.  2.3Validationofthenewmethodforestimatesoilhydraulicproperties The new method to estimate the soil hydraulic properties by inverse analysis of WCP was tested in three different porous media: 2 mm sieved loam soil, sand (particle size of 50-1000 m) and glass microspheres (particle size of 50-100m). The soil particle-size distribution of the different media was measured using the laser diffraction technique (COULTER LS230). One replication of soil particle-size distribution was performed per sampling site. Pre-treatment for the loam soil included the organic matter removing with hydrogen peroxide, soil shaking with a water dispersant solution and ultrasonic treatment. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 34 2.3.1 Measurement of soil hydraulic properties 2.3.1.1 Determination of water retention curve The water retention curve during a drainage process was measured with pressure TDR-cells (Moret-Fernandez et al., 2012). This consists of a 50-mm-long and 50-mm internal diameter stainless steel cylinder attached to a porous ceramic disc (bubbling pressure of 0.5 bar) and closed at the ends with two aluminum lids (Fig. 4). A 49-mm-long and 3-mm-diameter stainless steel rod, which longitudinally runs through the centre of the cylinder, constitutes the inner rod of a coaxial TDR probe. This rod is connected to the inner wire of a female BNC connector, which is glued onto the upper lid of the TDR pressure cell. The two elements, the stainless steel rod and the cylinder, form a cylindrical coaxial line of 49-mm length and 50-mm internal diameter. Two aluminum rings attached to several rubber joints hermetically close the lids of the TDR-cell against the stainless steel cylinder (Fig. 4). The TDR-cell is connected to a TDR cable tester (Campbell TDR100) by a 1.2-m long RG 58 coaxial cable of 50  nominal impedance, and the TDR signals are transferred to a computer that records and analyses the TDR waveforms using the software TDR-Lab V.1.0 (Moret-Fernández et al., 2010). The TDR volumetric water content are estimated using the Topp and Reynolds (1998) form.  Figure 4. Schematic diagram of the pressure head TDR-cell To set up the pressure TDR-cell the, 5 cm long stainless cylinders were filled with the tested porous media. Then, the stainless steel rod of the TDR-cell was inserted in the center of the stainless cylinder, and the top of the TDR-cell was hermetically closed by screwing the upper aluminum ring to the upper TDR-cell lid. In order to regulate the outlet water flow, a dry ceramic pressure plate was placed on the bottom lid of the TDR-cell. The stainless steel core plus the upper TDR-cell lid were attached to the ceramic disc. The system was finally hermetically closed by screwing the lower aluminum ring to the bottom lid of the TDR-cell (Fig. 5). Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 35 A first θ measurement was performed at air-dry soil conditions, which can be approached to a soil pressure head of 166 MPa (Munkholm and Kay, 2002). The soil samples were subsequently saturated by injecting distilled water through the base of the TDR-Cell, and the porous media were considered saturated when the water started to leave via the top of the pressure cell. Once the porous media was saturated, pressure steps were sequentially applied for 2-mm sieved loam soil at 0.5, 1.5, 3, 10, 50, 100, 500 and 1500 kPa. For the sand and the glass microspheres, additional pressure heads at 0.5, 5, 7.5 20, 40 and 70 kPa were supplied. Measurements of θTDR at the different pressure heads were done every 24 hours, except for the 500 had 1500 Kpa of pressure heads, where the samples were drained during 48 and 96 h. respectively. The parameters for the modeled unimodal and bimodal water retention curves (Eqs. 2 and3) were calculated from the experimental data using the SWRC fit (Seki, 2007) software. The dry bulk density of samples was also calculated from the volume and weight of the soil core after drying the samples at 105ºC during 24 h. Figure 5. Set of TDR-cell to measure the water retention curve In order to characterize the hysteresis phenomena on the sand media, an additional experiment was performed to measure the water retention curve during a wetting process. To this end, the TDR-cell, filled with dry sand, was connected by the top to an air pressure system, that inject air a constant pressure head. Next, distilled water was added by the bottom of the cell in that way that water raised by capillarity along the sand column while the pressure air was injected by the top. Once the water front arrived to the top of the soil column, the water content was measured by TDR. This process was repeated for an initially dry sand sample at 0.3, 0.5, 1 , 2, 3 ,4, 5, 6.5 and 500 kPa of injecting air. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 36 2.3.1.2 Determination of soil hydraulic conductivity and sorptivity The K and S parameters for a falling water infiltration process were measured using a minitension disc infiltrometer (Perroux and White, 1988) with 5 cm diameter. This instrument, made of Plexiglass, consists of a disc base covered by a membrane and a water supply reservoir (Fig. 6). The air inlet in the disc base was at 0.1-cm height from the soil surface. The cumulative infiltration curve, from which K and S are calculated, was measured from the water level drop of the water reservoir. This was measured with a 0.5 psi differential pressure transducer (PT) (Microswitch, Honeywell), that connected to a datalogger (CR1000, Campbell Scientist Inc.), was installed at the bottom of the water supply reservoir (Fig. 6). The infiltration experiment was made on an upside down TDR-cell head attached to a 5 cm long stainless cylinder filled with the selected porous material. Once the soil-filled cylinder was leveled, the mini-disc infiltrometer, filled with distilled water, was placed on the porous media core (Fig. 5). The level drop from the water at the reservoir tower was recorded at 5 s time interval until the water started to drain by the bottom of the cylinder. The K and S were numerically calculated (Latorre et al. 2013) by looking for the best fitting between the experimental and the theoretical quasi-analytical solution (Haverkamp et al., 1994) of the cumulative infiltration curve for disc infiltrometers. Figure 6. Mini tension disc infiltrometer Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 37 2.3.2 Experimental design to estimate soil hydraulic properties by inverse analysis of the water content profiles The WCPs under falling water infiltration were measured using the same TDR-cell head plus stainless steel core system described in section 2.3. For WCP measurements, the TDR-cell head was connected to a TDR100 (Campbell Scientific, USA) cable tester, which, using the TDR-Lab V.1.2. (Moret-Fernádnez, et al. 2011), transfers the TDR waveforms to a computer (Fig. 7). As described in section 2.3.1.2 the wetting process was performed from the top of the porous media core using the modified mini-infiltrometer (Fig 7). This system, allowed simultaneous measurements of WCP and cumulative infiltration curves. The TDR waveforms were recorded every 4 seconds, and the WCPs were calculated from inverse analysis of the TDR waveforms (see section 2.1). This data was subsequently treated with the Matlab-HYDRUS 1D interface to estimate the n, α, and Ks parameters. Finally, the n, α, and Ks parameters calculated with Matlab-HYDRUS 1D interface from the dynamics of the WCPs were compared to the water retention curve parameters for a draining process obtained with the TDR-cell experiment (section 2.3.1.1) and the saturated hydraulic conductivity calculated from disc-infiltrometer cumulative infiltration curve (section 2.3.1.2).  Figure 7. Experimental design for WCP measurement by TDR and cumulative infiltration during a falling water infiltration process Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 38 3.Resultsanddiscussion 3.1Experimentalresults 3.1.1 Particle size distribution of the porous media Different particle size distribution was observed in the different porous media (Fig. 8). The glass microspheres media, which showed the most homogeneous particle size distribution, contrasted with the sieved loam soil, with a gradient of particle size distribution. An intermediate state was observed in sand.    Figure 8. Particle size distribution of the three porous media    Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 39 3.1.2 Measurement and modeled water retention curve Experimental and modeled values of the water retention curve (WRC) for the 2-mm sieved loam soil, sand and glass microspheres are showed in Figure 9. The WRCs were fitted according to the unimodal (Eq. 2) and bimodal (Eq. 3) functions. In all cases, the WRCs show a better bimodal fitting, which indicates, that the employed porous media present a double-porosity architecture. Parameters for the WRC for a draining process fitted according to the uni- and bimodal functions (Eq. 2) (Table 2) show that the highest and lowest saturated water content corresponded to the loam soil and glass microspheres, respectively. This result agrees to that found in the literature, where coarse media presents lower total porosity (Hillel, 2003). The higher α value found in sand, which was similar to that obtained by Moret-Fernández et al. (2012), indicates this media allows retaining less water at near saturated conditions. These results contrast to those obtained in glass micro-spheres which needed applying high pressure heads to drop the water content at near saturated conditions. The lower n value in sand denotes that an abrupt water content drop occurs in a short pressure head interval. This value, however, was lower than that observed by Moret-Fernandez et al. (2012) in a similar experiment. This difference could be explained by a different packing of the sand samples. These results contrasts to that obtained for the loam soil and glass micro-spheres, where smoother water content changes regarding to the pressure head were observed. The n value for glass microspheres was similar to that obtained for a silty soil (Lipiec et al., 2007). The WRC measured in sand during a wetting process shows, compared to that for a draining process, higher α and n values. These differences should be attributed to the hysteresis phenomena (Hillel, 2003), defined as the difference in the relationship between the water content of the soil and the corresponding water potential obtained under wetting and drying process. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 40 Table 2. Water retention curve parameters modelled for the three porous media during a draining and wetting process using the unimodal function (Eq.2) with a residual water content equal to zero. θs (cm3 cm-3) α (kPa-1) n R2 Draining process Sieved loam soil 0.465 0.130 1.13 0.97 Sand 0.447 0.297 1.57 0.92 Glass microspheres 0.437 0.0015 1.55 0.98 Wetting process Sand 0.34 0.52 2.32 0.95 Table 3. Water retention curve parameters modelled for the three porous media during a draining and wetting process using a bimodal function (Eq.3) with residual water content equal to zero. θs (cm3 cm-3) α1 (kPa-1) n1 w α2 (kPa-1) n2 R2 Draining process Sieved loam soil 0.480 0.619 1.12 0.770 0.0002 2.75 0.999 Sand 0.423 0.164 6.69 0.694 0.0045 1.97 0.998 Glass microspheres 0.456 0.150 1.35 0.177 0.0008 1.90 0.999 Wetting process Sand 0.37 4.72 49.30 0.094 0.48 2.44 0.96 Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 47 Comparison between modeled and measured WCP needs that the initial time during wetting process to be exactly fixed, otherwise aberrant results may obtained. Except for the loam soil, a good fitting was observed between the WCP estimated by TDR and those modeled by the Matlab- HYDRUS 1D interface (Figure 15). The worse fitting between measured and modeled WCP in loam soil may be due to the unimodal instead of the bimodal water retention function used in the model.  3.2 MatlabHYDRUS1Doptimizationinterface 3.2.1 Estimation of hydraulic parameters (α, n, K) for porous media The α, n and K values obtained from the measured WCP by the optimization process calculated with the Matlab-HYDRUS 1D interface are showed in Table 5. Table 6 shows the relationship between the hydraulic parameters obtained from the TDR-cell experiment for unimodal and bimodal WRC functions and the corresponding values modeled by HYDRUS 1D for the three porous media. Overall, the relationship between modeled and measured hydraulic conductivity obtained for the three porous media was acceptable, near to one (Table 6). These results indicate that this procedure may be a feasible method to estimate K by inverse analysis of WCP. These results, however, contrast to those obtained for the WRC parameters, where modeled WRC shapes differ to that obtained from the experimental curves (Figure 16). Several reasons may explain these disagreements between the modelled and measured WRC parameters: 1. The HYDRUS 1D-Matlab uses a unimodal WRC model (Eq. 2), while the porous media fit better with a bimodal WRC function (Eq. 3) (Table 2 and 3). 2. Soil hysteresis phemomena (Hillel, 2003): while HYDRUS 1D calculates n and α parameters for a wetting process, the WRC obtained with TDR-pressure cells corresponded to a draining process. This result is visible in the sand experiment, in which the WRC for a wetting process is closer to the WRC obtained with HYDRUS 1D – Matlab interface (Figure 14b). As observed in Fig. 14a, the shape of the WRC for a wetting process is more similar to that calculated by HYDRUS, than the corresponding curve obtained during draining process. This hysteresis phenomenon may explain also the discrepancy between modelled and measured WRC on the sieved loam soil and glass microspheres (Figure 14a and c). To prevent this problem, hydraulic parameters calculated by HYDRUS 1D should be compared with the corresponding ones measured during a wetting process. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 48 3. Uncertainties of Eq. (52) to estimate the WCP from inverse analysis of the TDR waveforms. The testing of Eq. (52) for a same TDR waveform (data not shown) have revealed very small variations of the a and b coefficients can give different shapes of WCP. This problem could be solved by linking the HYDRUS 1D-Matlab interface with TDR physical model. This prevented using an intermediate function between both models, allowing the HYDRU 1D – Matlab interface to work with a single function of error that included the modelled TDR waveforms. Table 5. Hydraulic parameters modeled with Matlab-Hydrus-1D and root mean square error (RMSE) obtained for the three porous media.  α (kPa-1) n K (cm min-1) RMSE Sieved loam 0.01 3.29 0.001 2.25 Sand 1.47 5.50 1.402 0.06 Glass 11.00 2.95 3.281 0.59 Table 6. Relationships between hydraulic parameters obtained from experimental measurements fitted to the unimodal (Eq. 2) and bimodal (Eq. 3) functions and those modeled with the Matlab-HYDRUS 1D for the three porous media. αm/αexp n m/nexp K m/Kexp Unimodal Sieved loam 0.09 2.92 0.11 Sand draining 4.95 3.49 1.80 Glass - 1.90 2.23 Sand wetting 2.83 2.37 Bimodal αm/α1exp nm/n1exp Sieved loam 0.019 2.94 Sand draining 8.94 0.82 Glass 73.28 2.19 Sand Wetting 0.31 0.11 Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 49 (Kpa) 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 (m3m - 3) 0,0 0,1 0,2 0,3 0,4 0,5 0,6 a b (KPa) 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 (m3m - 3) 0,0 0,1 0,2 0,3 0,4 0,5 c (KPa) 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 m 3 m -3 ) 0,0 0,1 0,2 0,3 0,4 0,5 Figure 14. Modeled by HYDRUS 1D (blue continuous line) and water retention curves measured with the TDR-cell for a draining process on 2- mm sieved loam soil (a), sand (b) and glass micro-spheres (c). Triangles denote the water retention curve measured in sand during a wetting process.  Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 50 3.2.2 Sensitivity analysis of hydraulic parameters Figures 15, 16 and 17 show the results for sensitivity analysis of the function of error for α, n and K (see section 2.3.4). Functions of error for the sieved loam soil and sand, which present smoother shapes, allowed faster optimization techniques. For sieved loam soil and sand, the global minima found with the zoom optimization technique (Table 5), and the corresponding minimum values obtained by the sensitivity analysis (Fig. 15, 16 and 17) were quite similar. Although the sensitivity analysis for the α and K obtained for glass microspheres was acceptable, the function irregularities on the shape boundaries suggest that larger range to sweep the optimization process should be employed. The irregularities observed in the sensitivity analysis on n for the glass microspheres media suggest that the swept range is not adequate enough, or the initial conditions related to glass microspheres are not adequate.  Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 51 a  10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 Error 0 1 2 3 4 α min = 1.21E-03 Error = 2.21E-01 b  10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 Error 0 1 2 3 4 α min = 1,47 E-01 Error = 6,65E-02 c  10 -3 10 -2 10 -1 10 0 10 1 10 2 Error 0 1 2 3 4 α min = 1.1 Error = 8.47E-02 Figure 15. Sensitivity analysis for the  parameter calculated on sieved loam soil (a) sand (b) and glass microspheres (c) media.  Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 52 a n 1234567 Error 0,0 0,5 1,0 1,5 2,0 2,5 3,0 n = 3.31 Error = 2.25E-01 b n 1234567 Error 0,0 0,5 1,0 1,5 2,0 2,5 3,0 nmin = 5,38 – 5,5 Error = 6,65E-02 c  1234567 Error 0,0 0,5 1,0 1,5 2,0 2,5 3,0 n min = 2.95 Error = 0.0838 Figure 16. Sensitivity analysis for n parameter calculated on the sieved loam soil (a). sand (b) and glass microspheres (c) media.   Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 53 a K (cm/min) 0,0001 0,001 0,01 0,1 1 10 Error 0,0 0,5 1,0 1,5 2,0 2,5 3,0 K min = 1.27E-03 Error = 2.31E-01 b K (cm//min) 0,001 0,01 0,1 1 10 Error 0,0 0,5 1,0 1,5 2,0 2,5 3,0 K min = 1.4 cm/min Error = 6.65E-02 c K (cm/min) 0246810 Error 0,0 0,5 1,0 1,5 2,0 2,5 3,0 Kmin = 3.28 cm/min Error = 0.0917 Figure 17. Sensitivity analysis for the hydraulic conductivity (K) calculated on sieved loam soil (a). sand (b) and glass microspheres (c) Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 54  4.Conclusions This work presents a new TDR based methodology to estimate soil hydraulic parameters (α, n and K) by inverse analysis of water content profiles (WCP). The method is based on the RMSE analysis between the WCP estimated by the inverse analysis of TDR waveforms during a falling water infiltration process and the WCP simulated with HYDRUS 1D by sweeping a wide range of hydraulic parameters during a brute force optimization process. The method was tested on 2 mm sieved loam soil, sand and glass microspheres. The WRC parameters and the hydraulic conductivity obtained with this method were compared to those measured with conventional methods: TDR-pressure cell and disc infiltrometry technique. Although results show that the method allows satisfactory estimations of K, inaccurate n and α value were so far obtained. These results led to further efforts to improve the accuracy of the proposed method are needed. These new work should be addressed to: (i) employ bimodal instead of the unimodal WRC functions in the HYDRUS 1D simulations; (ii) link the HYDRUS 1D- Matlab interface to the physical model that TDR waveforms, and make a single function of error between modeled and measured TDR waveforms; (iii) employing faster optimization techniques. e.g. the Marquardt-Levenberg algorithm, Montecarlo algorithms or genetic algorithms. (iv) measuring WRC during a wetting process in order to take into account the hysteresis phenomena. On the hand, new efforts and experiments should be done to adapt this method to a capillary water absorption process, which would allow removing the gravity effect in the estimate of the soil hydraulic.       Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 55 References Angulo-Jaramillo, R, Vandervaere, J.P., Roulier, S, Thony, J.L., Gaudet, J.P., Vauclin, M. 2000. Field measurement of soil surface hydraulic properties by disc and ring infiltrometers. A review and recent developments. Soil Tillage Research 55: 1–29. Ankeny, M.D., Kaspar, T.C., Horton, R., 1988. Design for an automated tension infiltrometer. Soil Science Society of American Journal 52: 893–896. Ball J.A.R., 2002. Characteristic impedance of unbalanced TDR probes. IEEE Transactions on Instrumentation and Measurement. 51: 532-536. Brooks, R, Corey, A.T. 1964. Hydraulic properties of porous media. Hydrology Papers. Paper 3. Colorado State Univ., Fort Collins, CO, USA. Chong, E.K.P, Zak, S.H. 2013. An Introduction to optimization. John Wiley & sons. Third edition. ISBN: 978-1-118-51515-0. Connolly, R.D. 1998. Modelling effects of soil structure on the water balance of soil-crop systems: a review. Soil Tillage Research 48: 1–19. Dane J.H., Hopmans J.W. 2002. Water retention and storage. In Methods of Soil Analysis. Part. 4, Dane JH and Topp GC (editors). SSSA Book Series No. 5. Soil Science Society of America: Madison, WI. Durner, W. 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resources Research. 30: 211-223. Dexter, A.R. 2004. Soil physical quality: Part I. Theory, effects of soil texture, density, and organic matter, and effects on root growth. Geoderma 120: 227-239. Friel, R., Or, D. 1999. Frequency analysis of time-domain reflectometry (TDR) with application to dielectric spectroscopy of soil constituents. Geophysiscs 64: 707-718. Gavin H. 2011. The Levenberg-Marquardt method for nonlinear least squares curve fitting problems. Department of Civil and Environmental Engineering. Duke University. (http://people.duke.edu/~hpgavin/ce281/lm.pdf). Giese, K., Tiemann., R. 1975. Determination of the complex permittivity from thin-sample time domain reflectometry: improved analysis of the step waveform. Advances Molecular Relaxation Processes. 7: 45-49. Guérif, J., Richard, G., Dürr, C., Machet, J.M., Recous, S., Roger-Estrade, J. 2001. A review of tillage effects on crop residue management, seedbed conditions and seedling establishment. Soil Tillage Research 61: 13–32. Greco, R. 2006. Soil water content inverse profiling from single TDR waveforms. Journal of Hydrology 317: 325-339. Green, H.E., Cashman, J.D. 1986. End effect in open-circuited two wire transmission lines. IEEE Trans. Microwave Theory Techniques 34: 180-186. Model for soil hydraulic parameter estimation by TDR Carolina Peña Sancho 56 Haverkamp, R., Ross P.J., Smettem, K.R.J., Parlange, J.Y. 1994. Three dimensional analysis of infiltration from the disc infiltrometer. Part 2. Physically based infiltration equation. Water Resources Research 30: 2931-2935. Heimovaara, T.J. 1994. Frequency domain analysis of time domain reflectometry waveforms 1. Measurement of the complex dielectric permittivity of soils. Water Resources Research 30: 189–199. Heimovaara, T.J., Huisman, J.A., Vrugt J.A., Bouten, W. 2004. Obtaining the spatial distribution of water content along a TDR probe, Using the SCEM-UA bayesian inverse modelling scheme. Vadose Zone Journal 3: 1128–1145. Hillel D., 2003. Introduction to Environmental Soil Physics. Elsevier, Academic Press. Huebner C., Kupfer, K. 2007. Modelling of electromagnetic wave propagation along transmission lines in inhomogeneous media. Measurements Science. Technology. 18: 1147–1154. Kiefer, J. 1953. Sequential minimax search for a maximum. Proceedings of the American Mathematical Society. 4: 502-506. Kosugi, K. 1996: Lognormal distribution model for unsaturated soil hydraulic properties. Water Resources Research 32: 2697-2703. Latorre B., Motet-Fernández D., Caviedes D. 2013. Estimate of soil hydraulic properties from disc infiltrometer three-dimensional infiltration curve. Part 1. Theoretical analysis. Journal of Hydrology (submitted). Lera, F., Moret-Fernández, D., Vicente, J., Latorre, B., López, M.V. 2011. Comparison of different methods of TDR waveform analysis for soil water content and bulk electrical conductivity measurements using short-TDR probes. (to be submitted). Lin, C.P., Chung, C.-C., Huisman, J. A., Tang, S.H. 2008. Clarification and calibration of reflection coefficient for electrical conductivity measurement by time domain reflectometry. Soil Science Society of American Journal 72: 1033-1040. Lipiek J., Walczak R., Witkowska-Walczac B., Nosalewicz A., Slowinska-Jurkiewicz and Slawinski C. 2007. The effect of aggregate size on water retention and pore structure of two sily loam soils of different genesis. Soil and Tillage Research 97:239-246. Malicki, M.A., Plagge, R., Roth, C.H. 1996. Improving the calibration of dielectric TDR soil moisture determination taking into account the solid soil. Journal of Soil Science 47: 357-366. Meissner, T., Wentz, F.J. 2004. The complex dielectric constant of pure and sea water from microwave satellite observations. IEEE Transactions on Geoscienc and Remote Sensing 42: 1836-1849. Moret-Fernández, D., Vicente, J., Latorre, B., Lera, F., Castañeda, C., López, M.V., Herrero, J. 2012. TDR pressure cell for monitoring water content retention and bulk electrical conductivity curves in undisturbed soil samples. Hydrological Proccesses 26:246-254. Moret-Fernández, D., Vicente, J., Lera, F., Latorre, B., López, M.V., Blanco, N., González- Cebollada, C., Gracia, R., Salvador, M.J., Bielsa, A. Arrúe, J.L. 2010. TDR-Lab Version 1.2 User’s Guide. 2011.