Numerical study of the evaporation process and parameter estimation analysis of an evaporation experiment

Evaporation is an important process in soilatmosphere interaction. The determination of hydraulic properties is one of the crucial parts in the simulation of water transport in porous media. Schneider et al. (2006) developed a new evaporation method to improve the estimation of hydraulic properties in the dry range. In this study we used numerical simulations of the experiment to study the physical dynamics in more detail, to optimise the boundary conditions and to choose the optimal combination of measurements. The physical analysis exposed, in accordance to experimental findings in the literature, two different evaporation regimes: (i) a soil-atmosphere boundary layer dominated regime (regime I) close to saturation and (ii) a hydraulically dominated regime (regime II). During this second regime a drying front (interface between unsaturated and dry zone with very steep gradients) forms which penetrates deeper into the soil as time passes. The sensitivity analysis showed that the result is especially sensitive at the transition between the two regimes. By changing the boundary conditions it is possible to force the system to switch between the two regimes, e.g. from II back to I. Based on this findings a multistep experiment was developed. The response surfaces for all parameter combinations are flat and have a unique, localised minimum. Best parameter estimates are obtained if the evaporation flux and a potential measurement in 2 cm depth are used as target variables. Parameter estimation from simulated experiments with realistic measurement errors with a two-stage Monte-Carlo Levenberg-Marquardt procedure and manual rejection of obvious misfits lead to acceptable results for three different soil textures. Correspondence to: K. Schneider-Zapp (klaus.schneider@iup.uniheidelberg.de)


Introduction
Evaporation from porous media is a key process for soilatmosphere interaction, for example in the coupling with climate or the forcing of lower soil layers, as well as for many industrial and engineering applications.Many investigations are reported in the literature which assess the evaporation process.Evaporation from an initially saturated porous medium typically begins with a relatively high drying rate determined primarily by the external forcing.This phase continues as long as the medium can sustain the evaporative flow.Then it changes to a stage with falling drying rates (Sherwood, 1930;Scherer, 1990;Shokri et al., 2008).Extensive work has been put into pore-scale modelling of the drying process; a review is given by Prat (2002).The pore-scale analysis is valuable for the understanding of the detailed pore-scale processes.However these models typically cannot be directly applied for macroscopic problems, since the actual geometry of the medium is normally not known.On the field scale, many energy-balance based, semi-empirical, or empirical models exist (Foken, 2003).For modelling on the REV scale, in-between the field-and pore-scale, Schneider et al. (2006) used a diffusive boundary layer approach coupled with a Richards' pore space model.It is reasonably simple but still provides a sufficient macroscopic description.One aim of this investigation was to further examine the physical implications of that model.
Measurements of the hydraulic properties of soils in the dry range are hard to realise.While direct measurements of hydraulic properties are generally difficult, Multistep-Outflow experiments are limited to matric potentials ψ m >−100 kPa since the liquid phase pressure must be larger than the vapour pressure of water.Practical limitations like the permeability of the phase separator at the lower boundary are more strict and typically lead to ψ m >−20 kPa.
Published by Copernicus Publications on behalf of the European Geosciences Union.Similar restrictions apply to traditional evaporation experiments, where a saturated soil sample is placed on a balance and exposed to free air while the matric potential in several depths is measured by tensiometers.If the potential falls below the air-entry value of the tensiometer or below the vapour pressure of water, whichever is higher, then water is released from the tensiometer into the soil.This leads to a disturbance of the measurements that may be quite dramatic.The measurement range is further limited by the technical challenges to measure with tensiometers the very small potential gradients in regions where the hydraulic conductivity is still high, or to assess the little weight change caused by small evaporation fluxes in the dry range with a balance.Schneider et al. (2006) presented a novel evaporation experiment which yields reliable data even in the dry and very dry range by measuring the flux at the upper boundary with an infrared absorption gas analyser.An inverse model for the estimation of hydraulic parameters from the evaporation flux was developed.The analysis of an evaporation experiment with an undisturbed soil sample yielded reasonable results.
The objective of this study is to analyse the properties of this novel evaporation experiment in more detail by conducting virtual experiments and perform parameter estimation on this synthetic data.Specifically, we (i) use the model to study the physical processes during the experiment, (ii) study the sensitivity of the measured quantities to parameter changes to optimise the boundary conditions of the experiment, (iii) explore if adding more observables into the inversion process obtains significantly more information about the system, and (iv) analyse the identifiability and uniqueness of the solution.
Multi-dimensional non-linear optimisation problems often have more than one minimum and the minima are often not well-localised leading to ambiguous or contradictory solutions.Thus it is important to preclude such a behaviour with a detailed analysis.

Setup of the novel evaporation experiment
The experimental setup of the novel evaporation experiment is described in detail in Schneider et al. (2006).Therefore only a brief description is given here.
The soil sample is contained in a PVC cylinder of 10 cm height.The bottom of the column is closed, the top of the soil column is closed by a gas-tight head space (evaporation chamber) (Fig. 1).A constant flow of air is established through the head space to remove the water vapour.The water vapour partial pressure p w and temperature T of the incoming air are controlled and thereby the boundary condition at the upper boundary is set.The air in the head space is turbulently mixed, leading to a uniform potential in the chamber.The water flux is quantified by the difference of water vapour content before and after the evaporation chamber and the prescribed air flow through the head space.

Numerical model
We model the soil column as a uniform one-dimensional medium and assume that its soil water characteristic may be described by the van Genuchten parametrisation and its hydraulic conductivity function by the corresponding Mualem parametrisation The transport model is developed in detail in Schneider et al. (2006).Assuming local thermodynamic equilibrium, the molar water vapour content ν w g [mol m −3 ] is given by the Kelvin equation (Rawlins and Campbell, 1986) where V w m = 1.804 × 10 −5 m 3 mol −1 is the molar volume of liquid water, R = 8.3145 J mol −1 K −1 the universal gas constant, and p w s (T ) [Pa] the saturation partial pressure of water vapour over pure liquid water at temperature T [K].It can be described with Magnus' formula (Murray, 1967) These relations can also be used to calculate the equivalent matric potential from a given water vapour concentration.
The equivalent flux j w g [m s −1 ] of liquid water transported by diffusion of water vapour is given by where D w g [m 2 s −1 ] is the diffusion coefficient of water vapour in air.By using the chain rule for ∇ν w g , neglecting the temperature dependent part of water vapour diffusion in soil, approximating water vapour by an ideal gas, and describing water vapour diffusion by D w g =θ 4/3 s D w g,atm (Jin and Jury, 1996), where D w g,atm is the diffusion coefficient for water vapour in free air, we can include vapour transport in Richards' equation as an effective conductivity: A crucial step is the representation of the upper boundary.The transition between soil and atmosphere will affect the potential, resulting in in a potential in the upper soil which is higher than the effective potential of the atmosphere.The air in the head space is turbulently mixed, leading to a uniform potential throughout the chamber.The eddies cannot penetrate the boundary, thus eddy size will decrease when approaching the surface.In a thin layer around the surface, diffusive transport will dominate.Since turbulent mixing is much more efficient than diffusive transport, the vapour flow is controlled by a thin diffusive layer around the soilatmosphere boundary.Because the soil is rigid, the boundary layer thickness is not likely to change.We assume that the time scale of diffusion across this layer is much smaller than the time scale on which the boundary condition changes.This appears reasonable since the time scale of diffusion, given by r 2 b /[2D w g,atm ], is some 0.1 s for a layer thickness of 2 mm.The vapour flux across such a layer is given by where p w exp [Pa] is the partial pressure of water vapour in the atmosphere, T [K] the temperature in the well-mixed head space above the soil column, and r b [m] the effective thickness of the boundary layer.We comment that, by definition, the processes in this layer are not resolved well.In particular its physical location is not defined, i.e., the fraction of the layer that is within the soil column, and the porosity of the respective parts of the soil.However, for the type of thin layer we consider here all these complicating factors only enter as a constant of proportionality.Hence, we have made r b a fitting parameter that absorbs all these factors.Obviously, its value then cannot be interpreted physically any more.
Hydraulic parameters were estimated from the measurements using inverse modelling.As it is impossible to estimate both θ s and θ r by an evaporation experiment alone, we set θ r to zero for all experiments and just fitted θ s , which then corresponds to a kind of "available water content".Further fitted parameters are the Mualem-van Genuchten parameters α and n, the saturated hydraulic conductivity K s , and the effective resistance r b of the boundary layer.The value of τ was fixed at 0.5 as suggested by Mualem (1976).
We used the numerical forward model together with the Levenberg-Marquardt algorithm (Marquardt, 1963), where the residuum is calculated by the squared sum of normalised deviations: where y measured i denotes the ith measured data value, y model i the simulated measurement, σ i the standard deviation of the ith measurement and N the number of measurement points.The forward model integrates Richards' equation using a cell-centred finite-volume scheme with full-up-winding in space and an implicit Euler scheme in time.Linearisation of the nonlinear equations is done by an inexact Newton method with line search.The linear equations are solved with a direct solver.For the time solver the time step is adapted automatically.A no-flux condition was used for the lower boundary.At the upper boundary the evaporation was calculated by Eq. ( 8).We did not simulate energy loss due to the latent heat of evaporation and the heat transfer in the soil sample, assuming that the heat exchange between the sample and its environment is fast enough to compensate for the heat loss by evaporation.As shown in Schneider et al. (2006) this assumption is violated in the initial phase of the experiment with a sandy loam.While we do not expect principle differences in the system behaviour, the consequences of this simplification still have to be studied in the future.
The sensitivities required by the Levenberg-Marquardt algorithm were derived by external numerical differentiation.

Numerical experiments and analysis
Virtual experiments have been conducted where the true parameters are known and therefore the performance of the inversion process as well as the parameter space can easily be analysed.In our study, we used three parameter sets: a sand, a silt and a sandy loam.The corresponding parameters are given in Table 1.All experiments were simulated isothermally with T =293 K. Boundary conditions Table 1.Parameters used for the synthetic data sets: The van Genuchten parameters α and n, the saturated hydraulic conductivity K s , residual water content θ r , the saturated water content θ s , and the resistance of the boundary layer r b .parameter sand sandy loam silt which were finally used are given in Table 2.All experiments were started with a fully saturated soil sample.
To test if additional measurements can improve the quality of the parameter estimation besides the evaporation flux j w two additional (virtual) measurements are considered in this paper: (i) the matric potential ψ m measured by a tensiometer, and (ii) the liquid water content θ measured by the dielectric (compound) permittivity c [-].Both probes are assumed to be installed 2 cm below the surface.The influence of the installation depth is analysed below.The measurement uncertainty assumed for each of the virtual devices is given in Table 3.Note that, in contrast to real measurements, the virtual instruments provided point measurements.We calculated the square root of the permittivity depending on the water content according to the formula (Roth et al., 1990): where w is the permittivity of water, s of the soil matrix, and a of air, respectively, and β = 0.5.The saturated water content θ s was used as porosity and s was assumed to be known.
If the potential falls below the air-entry value of a tensiometer or below the vapour pressure of water, whichever is higher, the tensiometer releases water to the sample.To prevent this disturbance, the tensiometer is removed at −30 kPa.To analyse the impact of this removal we also conducted simulations were the tensiometer was removed at −70 kPa, and simulations with a (hypothetical) unlimited measurement range.
Measurements were made after logarithmically growing time intervals t i =300 s+2000 s×log(i), starting again with i=1 after each change of the boundary condition.
As large gradients are encountered in the simulated soil, especially at the drying front, a fine grid and small time steps are needed to avoid numeric noise.On the other hand, to keep the runtime reasonable, the spatial grid should be as coarse as possible.A grid convergence study was conducted showing that a reasonable grid convergence was obtained with a 1000 point non-regular grid with exponentially decreasing cell heights towards the soil surface.The uppermost cell had a height of 10 −9 m.Of course, this exceedingly small size is not related to the real physics at that scale.However, one has to bear in mind that the necessary grid resolution can also depend on the hydraulic parameters used in the simulation.The time step was adopted automatically by the model.
Forward simulations of the evaporation experiment were used to study the physical dynamics of the system.To obtain the maximum amount of information about the unknown parameters to be optimised, a sensitivity analysis was performed for the experiment.Relative sensitivity coefficients were calculated according to where m i denotes measured quantity i (e.g. the water flux at the upper boundary j w ) and p j is the value of the j th parameter.s i is a dimensionless quantity normalised by the measured quantity and parameter value which allows to compare the sensitivities of different measured quantities and different parameters, and is (except for numeric noise caused by the numeric differentiation) independent of the step size p j .It is a measure how much the quantity i changes when changing parameter j .When s ij is positive, quantity i increases when increasing parameter j , when it is negative, quantity i decreases.The larger the absolute value of s ij , the higher the change of quantity i for a given parameter change.When encountering a zero crossing of s ij , quantity i changes in opposite directions before and after the crossing, while it is not affected by parameter changes directly at the zero point.The results of the sensitivity analysis were used to optimise the boundary conditions of the experiment.
To check whether the data measured in the experiment is sufficient to identify a unique set of soil hydraulic parameters, response surfaces were calculated for all scenarios (similar to Toorman et al. (1992) for the onestep outflow experiment and Šimůnek et al. (1998) for traditional evaporation experiments).Two parameters were varied independently while all other parameters were kept at their true values.The χ 2 surface, defined by Eq. ( 9), is then displayed in contour plots.The true value of two parameters j and k were multiplied by factors γ j and γ k , p j,test = γ j p j,true .The factors γ were varied from 0.1 to 2 in steps of 0.1 and from 2.2 to 3.8 in steps of 0.4, respectively.The χ 2 value corresponding to a certain parameter combination was then colour-coded in a (γ j ,γ k ) plot.If one localised minimum with spherical isolines exists in the contour plot this indicates a well posed problem.Valleys in the χ 2 surface indicate correlations between parameters, since χ 2 stays relatively constant although two parameters were changed.This makes it more difficult to identify a unique set of parameters.Multiple minima would indicate a non-unique solution.While this approach only shows a subset of the true five-dimensional parameter space along parameter planes and new features might occur in the intermediate space, it is nevertheless a good indicator whether one unique and identifiable minimum exists.Response surfaces were also used to investigate if adding more observables yields substantially more information for the inversion process.
Finally, it was checked how well the inverse model converges to the real parameters given the measurement error and the cross correlation between the model parameters.The forward model was used to generate synthetic data for the best combination of observations as determined from the response surfaces.Random noise normally distributed with a standard deviation of σ i was added to the data, where σ i is the measurement uncertainty in the evaporation experiment (Table 3) as determined in Schneider et al. (2006).Five different data sets were generated for each scenario to account for the random influence of measurement noise.From each of this data sets, the parameters α, n, K s and θ s of the van Genuchten/Mualem model and the resistance of the boundary layer r b were estimated with a variety of initial conditions.10 parameter sets were randomly created in the range reasonable for the soil under examination.A logarithmic random distribution between the upper and the lower limit was used.These sets were used as start parameters for the gradient based inversion process resulting in 50 sets of estimated parameters for each soil.This approach combines a Monte-Carlo method with the Levenberg-Marquardt minimisation.We call this a Monte-Carlo Levenberg-Marquardt method, MCLM.The inverse solutions were then compared with the real parameters and the resulting hydraulic functions with the true functions.A deviation coefficient was defined according to where p inv denotes the inverted and p true the true parameter.

Results
As the results are quite similar for the three different soil types, we discuss only the results for the silt in detail.The results of the inverse modelling are given for all three soil types.

Physics of the process
The most simple scenario for an evaporation experiment which is also used in classical evaporation experiments is a onestep experiment as shown in Fig. 2.After saturation the sample is exposed to a constant vapour pressure at the upper boundary resulting in a progressive drying of the sample.
In this scenario two different regimes are distinguishable: Regime I, where the outflow is limited by the resistance of the boundary layer r b , and regime II, where it is limited by the soil hydraulic properties.These two regimes are in accordance with experimental findings in the literature.Regime I leads to a nearly constant value of j w which mainly depends on r b .Close to saturation changes of potential only lead to small changes in the water vapour pressure at the soil surface and thus of the resulting flux.Therefore, in this regime hydraulic properties cannot be determined with only the outflux as measured quantity.A sketch of the potential profile (1) at an early time in regime I, (3) directly before the transition to regime II, (4) a short time later after the transition, (6) at the end of the experiment.Notice the non-linear scaling of the depth axis.near the surface for three times is shown in Fig. 3.The major part of the potential drop is caused by the resistance of the boundary layer.Due to the much higher conductivity in the soil, water is delivered to the evaporating surface with a minimal gradient, which can also be seen in the simulation results (Fig. 4a).Hence, the hydraulic properties have only a minimal influence on the flux.However, the potential of the sample changes with time due to the successive drying of the soil.Measurements of potential could therefore give information about the hydraulic properties of the sample.For t 20 h, the potential in the soil is above −10 kPa and therefore can be measured easily with a tensiometer.However, due to the small deviations from the linear decrease (which would be expected for a hydraulic conductivity which is constant over the whole sample), one would need a very high accuracy to obtain information about the hydraulic conductivity.As the relative permeability of the silt decreases only slowly with decreasing water potential, a small potential gradient is sufficient to sustain the water flux to the soil surface and the sample dries more or less uniform over the whole whole 10 cm high sample.
With continuing evaporation the potential and the water content decrease (Fig. 4), most rapidly near the surface.Eventually the conductivity of the soil becomes limiting.Shokri et al. (2009) found at this point in their experiments the formation of a dry surface layer, which is in excellent agreement with our simulations.The system enters regime II and j w starts to decrease rapidly.This transition is quite abrupt because (i) the function K(θ ) is very steep in the relevant range and (ii) the effective hydraulic conductivity is dominated by the dry low-conductive layers where vapour diffusion is the relevant process.Therefore the hydraulic properties seen in this regime are the properties of the dry region.During this transition a drying front (interface between unsaturated and dry zone with very steep gradients) forms at the surface and then moves into the soil (Fig. 4).Since the low-conductive layer (the distance inside the soil to be passed by vapour diffusion) grows with continuing evaporation, the flux continues to decrease.
Regime I occurs only if the saturated hydraulic resistance (the reciprocal of the saturated hydraulic conductivity), which is the lowest possible resistance inside the soil, is lower than the resistance of the boundary layer.This is illustrated by reducing K s by a factor of 20, resulting in K s = 0.05 cm/h (Fig. 5, blue and dashed cyan curves).K s is now lower than the flux which can be evaporated through the soil-atmosphere boundary layer.Thus, the sample directly enters regime II.

Sensitivity analysis
The relative sensitivity coefficients according to Eq. ( 11) were calculated for each measurement type for all hydraulic parameters (Fig. 6).
As the system is in regime I at the start of the experiment, the water flux is limited by the resistance of the soilatmosphere boundary layer alone.At early times, the outflux j w is therefore most sensitive to r b while the sensitivity to all other parameters is very small.
When the topmost layer of the soil has dried out, the soil hydraulic properties, in particular the hydraulic conductivity, become limiting for the evaporation rate.The evaporation flux is most sensitive to all parameters exactly at the bend point where the system enters regime II and the outflow starts to decrease after the plateau.The sensitivity on r b decreases continuously because the drier the sample, the less important is the resistance of the boundary layer.As θ s scales the amount of available water when θ r is held constant, the sensitivity of θ s stays more or less constant after a quick decay.For all other parameters the sensitivity decreases after the maximum and after a zero-crossing eventually increases again with opposite sign.The zero-crossing can be explained by mass conservation.As the total water content of the sample is constant, a higher evaporation at earlier times has to be compensated by lower evaporation towards the end of the experiment and vice versa.
In contrast to the evaporation flux, the potential ψ m is at the beginning of the experiment most sensitive to α and n which control the shape of the soil water retention curve.K s is the only parameter for which the sensitivity is nearly zero during regime I.The sensitivity on r b and θ s are less important at the beginning, but increase during regime I and reach a maximum at the transition to regime II as well as the sensitivity on n.
Since r b determines the speed of drainage in regime I, it is clear that the potential, which is connected to the water content by the water characteristic, is also dependent on r b .This effect will become more pronounced as time passes because the longer a different outflux caused by a different r b is re- Fig. 5. Water flux at the upper boundary j w and potential in 2 cm depth for the "onestep" experiment with 20 times lower K s , as well as for the twostep experiment with 20 times higher α.The former stays virtually no time in regime I due to the very low conductivity.The latter shows a more distinct switch-back to regime I but structurally it is identical to the original soil.Notice that p w jumps from equilibrium to 1 kPa and 0.25 kPa, respectively, at t=0. tained, the higher are also the differences in the potential.θ s determines the amount of available water.With a constant evaporation rate, the more water is available, the less is the relative change of water content and therefore the change of potential when all other parameters are kept constant.Thus, the behaviour of θ s is analogous to the one of r b .
The maxima are much less pronounced than for the evaporation flux.For all parameters the sensitivity is more or less constant or increases slowly in regime II and reaches a large peak at t≈370 h which is caused by the passing of the drying front at the tensiometer position.This peak is discussed in more detail in Sect.3.2.2.A zero-crossing only exists for n as only n influences the shape of the soil water retention curve.The generally higher sensitivity in the dry range is in accordance with the result of Šimůnek et al. (1998) for traditional evaporation experiments.As they pointed out, the water characteristic becomes steeper for more negative potentials and thus parameter changes have more influence at lower potentials.
The water content is less sensitive to parameter changes than the evaporation flux and the potential.During regime I the only sensitive parameters are θ s and r b , which both increase with time.The maximum is again at the transition to regime II.For all other parameters there is no pronounced maximum at the transition point, but all sensitivity curves show a sensitivity maximum at the passing of the drying front.The water content is most sensitive to n and K s during the early stage of regime II and to the available water θ s towards the end of the experiment.A detail view of the first 100 h for the potential in 2 cm is shown in (d).Notice that p w jumps from equilibrium to 1 kPa at t=0.
The sensitivity to changes of the saturated hydraulic conductivity is rather low for all types of measurements and reaches significant values only for the evaporation at the transition point and at the passing of the drying front.

Physics of the process
As the transition from regime I to regime II contains much information, one would suggest that multistep experiments can drastically improve the sensitivity if it is possible to reproduce the switch from regime I to regime II with boundary condition steps.To switch from regime II back to I, either the conductivity in the upper soil or the resistance of the boundary layer must be increased, or the potential drop on the boundary layer decreased.As the resistance of the boundary layer r b is constant, this cannot be achieved by lowering the boundary potential.Lowering the boundary potential speeds up the drainage of the sample but does not lead to new features.
If the sample is already in regime II and the water vapour pressure at the surface is increased, a second plateau and a second drop of the flux can be seen if the pressure jump and the time between the steps are chosen adequately.This is illustrated with a twostep experiment, Fig. 5 (red and dashed magenta curves).To make the second step more pronounced, a 20 times higher α was used for this simulation.When the vapour pressure at the boundary is increased, the potential drop over the boundary layer and thus the water flux decreases, the boundary layer becomes limiting again.The flow inside the soil is now higher than the evaporation flux.This leads to an increase of the water content and thus the Hydrol.Earth Syst.Sci., 14, 765-781, 2010 www.hydrol-earth-syst-sci.net/14/765/2010/ potential at the soil surface, resulting in a larger potential drop on the boundary layer and therefore an again higher evaporation flux (see Fig. 7).When this adaptation stage is finished, the increase of water content in the upper soil and therefore the increase of the evaporation rate ends, the delivery from below and the evaporation flux are equal again.If the change in the boundary condition was large enough, the resulting evaporation flux at this point is low enough to be sustained by the soil for a longer time span and regime I is reached again, else the system stays in regime II.This depends on the relation between the new potential drop on the boundary layer and the hydraulic conductivity in the soil (soil water state).
The same effect also occurs with the normal value of α, but it is harder to see (after the first boundary condition jump at t = 62.5 h in Fig. 8 red line).The higher α results in a less negative potential in the soil before the switch and a relaxation to a higher water content after the switch.Therefore it takes longer until the conductivity drops low enough to reach regime II again.
We acknowledge that any change in the direction of flow leads to hysteresis, which was not considered in the simulation.However, hysteresis is most pronounced in coarse textured porous media and in the wet range of experiments.In our numerical multistep experiments, the switch occurs only in the dry range where water vapour is the only relevant transport mechanism.Thus hysteresis should be negligible.
The twostep experiment has the disadvantage that the potential range covered is too small during a reasonable measuring time.This can be compensated by applying a third step after the experiment has entered the hydraulically dominated regime again to speed up drainage.This results in a threestep experiment (Fig. 8, red and dashed magenta curves) which has identical general features as the twostep experiment but has a much larger potential range in 2 cm depth during the same measuring time.
The same boundary conditions may not produce a second regime switch with any type of soil.If the boundary conditions are changed too early, the plateau just changes its level as the potential difference changes.In this case the threestep experiment does not enhance the estimation of the hydraulic properties.However, it still gives more information about the resistance of the boundary layer.Prior knowledge about the soil hydraulic properties is required to choose the optimal boundary condition steps in a multistep experiment.This knowledge could be obtained by first performing a onestep experiment and then using this information to design a multistep experiment.However, a major disadvantage of this scheme is the long time required to conduct two experiments.
To study the importance of water vapour flow inside the soil compared to the flow of liquid water, we also performed a simulation where the effective conductivity contributed by water vapour flow was disabled (Fig. 8).Water vapour flow inside the soil is especially important at later times, when the soil becomes very dry after the third step of the boundary condition at t=196.3 h.Without the water vapour transport the hydraulic conductivity is already too low to get an increase of the evaporation flux when the vapour pressure at the boundary is reduced.The sample is effectively sealed by a very dry layer at the sample surface with very low conductivity, which prevents the further drying of deeper regions.A second feature not present in the simulation without vapour transport in the soil is the "undershoot" of the evaporation at the transition back to boundary layer dominated regime at t=62.5 h.With vapour transport, the evaporation before the switch is higher and thus the potential at the surface is lower.This results in a more pronounced drop of the evaporation and a longer time till the dynamic equilibrium in the soil is reached again.

Sensitivity analysis
While there is no big change in the relative sensitivities for the potential and the water content 2 cm below the soil surface compared to the onestep experiment, the sensitivities of j w increase substantially for experiments that re-enter regime I (Fig. 9).Multistep experiments which do not reenter regime I show no strong effect in the sensitivities (data not shown here).Thus, multistep experiments are a good tool to increase the sensitivity.Considering the larger potential range which is covered by the threestep experiment, it is also considered superior to the corresponding twostep experiment.
To determine the optimal position of a tensiometer or permittivity probe, profiles of relative sensitivity have been analysed.Figure 10 (top) shows profiles of the relative sensitivity of the matric potential to changes in α.Generally, the sensitivity is lower near the sample surface especially at later times.After the onset of regime II, a large sensitivity peak appears, which moves downward with time.At t=44 h the sensitivity drop below the peak even leads to a zero-crossing of the sensitivity.After the increase of water vapour pressure at the upper boundary and the transition back to regime I the sensitivity peak vanishes for a short time and reappears after the fall-back to regime II.This is the result of a rewetting of the top soil, which is caused by the fluxes from below which are higher than the (reduced) evaporation rate.
The peaks are located at the drying front (as can be seen for t=176.9h from a comparison with Fig. 11) where the gradients are particularly high and thus small deviations in the parameters lead to large changes in the solution.As changes of the parameters also affect the position of the drying front, small parameter changes generate huge potential differences in its proximity.Figure 11 illustrates the change in the profiles of matric potential and water content at t=176.9 h if α is increased by 10%.A zero-crossing of the sensitivity occurs if the parameter change leads not only to a different position, but also to a change in the steepness of the drying front.This is in accordance to Romano and Santini (1999) who reported that sensitivities of the matric potential in traditional evaporation experiments in the uppermost part of the soil show increasing curvatures and a drop to zero.They noticed that this especially happens at larger times when " h changes its sign in close proximity to the evaporating surface", where h is the change of matric head.A change in the sign of h corresponds to a change of the sign of the corresponding sensitivity coefficient as it was found in our study.
The sensitivity of the water content to changes in α is similar (Fig. 10, bottom).However, there is always a relevant  sensitivity at the soil surface and a zero-crossing which is located above the peak at the drying front for later times.
In principle the profiles of the relative sensitivity are similar for all parameters.As shown in Fig. 11 the profiles of the sensitivity of the matrix potential at t=176.9 h for different parameters mainly differ in the size and the sign of the peak.The sensitivity to changes of n has a very pronounced zero-crossing as n always influences the steepness of the drying front.The sensitivity peaks are most pronounced for the parameters θ s and r b as these parameters have the strongest influence on the propagation speed of the drying front, by determining the speed of drainage in regime I and thus the starting time of the drying front movement.
For the permittivity probe a position nearer to the surface than the 2 cm we used would be advantageous as the sensitivity of the water content is always high there, but this is hard to realise experimentally.For the tensiometer a depth of 2 cm is quite fine, as the sensitivity at this depth is high except for very late times after the passing of the drying front.When  Relative sensitivity s(z,t) of the potential ψ m on all variables at t=176.9 h (top).The potential and water content profor the undisturbed parameters as well as with α changed by 10 % (bottom) illustrate the source of the huge sensitivity peaks at z≈0.15 cm in the plot above, the location is marked with a dotted vertical line.Note the non-linear scaling of the depth axis.
the drying front passes, the potential drops so low (−10 4 kPa -this corresponds to −1 km water column -or less, see Fig. 11, bottom) that it is outside the measurement range of traditional tensiometers.Thus, for a real measurement tensiometers have to be removed before the drying front passes to avoid a leakage of water into the soil and therefore the sensitivity peaks of potential cannot be utilised with traditional tensiometers.
Especially the sensitivity peaks of the water content measurements open up the possibility to "scan" different sample layers during one experiment, as the sensitivity is focused on a very small height interval and penetrates with time.If the soil has different layers, the sensitivity would penetrate through these layers with the drying front.However, measurements of the whole water content profile would be necessary to exploit this.If only point measurements are available, the gathered information is not sufficient to distinguish influences at different depths, i.e. if changes in the drying front propagation are caused by another layer above or by different parameters of the same layer.The measurement of water content profiles could be done e.g. by X-ray, neutron or gamma adsorption.K. Schneider-Zapp et al.: Numerical analysis of evaporation experiment Fig. 12. Response surfaces for the threestep experiment with j w and ψ m as target variables.The light regions with deviation coefficients larger than 2 were calculated with a smaller resolution of 0.4.Non-matching structures between the two regions are caused by the different resolution.Note that a χ 2 above 10 7 is displayed in white, points where the model did not finish after 120 000 time steps is marked grey.

Response surfaces
Figure 12 shows the response surfaces of the threestep experiment with j w and ψ m as target variables.The tensiometer was removed at −30 kPa.Generally, there is a single global minimum and it is relatively well-defined.Only the combinations (α,K s ), (n,K s ), (K s ,θ s ), and (K s ,r b ) where K s is involved have small valleys and the slope to the absolute minimum in the direction of the valleys is low.This is a conse-quence of the low sensitivity of the measurements to changes in K s .As n=2 and the van Genuchten/Mualem parametrisation does yield non-physical hydraulic conductivity functions for n<2 (Ippisch et al., 2006), the part of the response surface with n<2 (i.e.n/n 0 < 1 in the contour plots) must be regarded with care.For the other target variable combinations and boundary conditions the response surfaces look similar and thus are not all shown here.To identify the essential measurements for a good estimation of the parameters the response surfaces of the threestep experiment for K s and n for combinations of the three measurement types are analysed (Fig. 13) as K s is especially hard to estimate due to its low overall sensitivity.If only the evaporation flux j w is used the residual does not have a well defined minimum but more a banana like shaped extended region.The minimum in the direction of n is better defined if only the matric potential in 2 cm depth is used but there is still a very long valley in the direction of K s .This valley becomes much shorter if a combination of j w and ψ m is used.The addition of permittivity (i.e. water content) measurements improves the situation only in the combination c +ψ m , in all cases involving the evaporation flux the changes are very small.However, it should be noted that accurate permittivity measurements (e.g. with TDR probes) in the dry range are not feasible.The reasons are: (i) Since the traveltime error is constant, the relative error increases with decreasing c .(ii) The compound permittivity c is a function of the soil matrix permittivity s , of the porosity φ, and the actual geometry.For low water contents, the uncertainty of the measurement diverges since the permittivity contribution of the remaining water becomes equal to or even lower than s , φ is not known accurately, and the geometry is unknown.(iii) For thin films of water as found in the dry region, w is different from the one of bulk water.Additionally, the measurement volume which was neglected in the simulations will smear out gradients, and results are expected to be worse than in the ideal case of a point measurement.
As j w is the derivative of the total water content, it is reasonable that adding a water content measurement does only give slightly more information.In contrast, ψ m is an independent observable and therefore gives more information about the soil water retention curve.However the fundamental difficulties with tensiometers must be regarded.The information is only provided in a small potential window and great care must be taken that the tensiometer is removed before the potential in the soil drops below its air entry point, or the water vapour pressure, whatever is higher.Therefore it was also investigated whether a tensiometer which can measure up to ψ m =−70 kPa or even without a limit makes a con- Comparison between the accepted fits with the largest residuum and the rejected fits with the lowest residuum, for sand (a) (residuum 4900 (accepted) vs. 9700 (rejected)), loamy sand (b) (residuum 28 000 vs. 3300), and silt (c) (residuum 4600 vs. 7300), respectively.The reason for rejecting the curves were for (a) the deviation of the tensiometer, for (b) the first plateau was not represented, and for (c) the first peak was not represented, respectively.Note that for (b), although the potential of the accepted fit matches worse than in the accepted one, its deviation is just at the sort-out limit and limit cases were still accepted, while the plateau in the rejected fit is definitely not represented at all.siderable difference.The results on the response surface are illustrated in Fig. 14.The lower the tear-off of the tensiometer, the better the data.Going from ψ m =−30 kPa to −70 kPa gives a significant enhancement.For the hypothetical case of unlimited potential measurement, the minimum is so localised that χ 2 is above the upper colour scale limit for all values but the minimum itself (data not shown).Apparently, tensiometers with a much wider range of measurements like the ones presented by Bakker et al. (2007) would be most helpful for this type of experiment.

Convergence study
The convergence study revealed the existence of local minima which are not visible in the response surfaces, because they are not located in two-dimensional sub-planes of the five-dimensional parameter space.When initially running the inverse fits with potential and outflux data, the convergence was poor.The reason was that when the outflux peaks did not fit initially and the algorithm slightly modified the so-lution vector to improve it, the potential in the new solution fitted worse.Because the potential has a relatively low standard deviation, the residuum was not improved and therefore the inverse fit could not improve the solution.To resolve that problem, the inversion was first run with the outflux data only to obtain a reasonable start parameter set.After convergence, the potential data was added and the inversion was restarted.
Using that modified approach there were still a few nonconverged fits.However these could be identified clearly because the system response was apparently not fitting the data.These fits were taken out manually.The criteria for sorting out were (1) the outflux peaks were not represented, (2) the potential systematically deviated by more than 10σ Hydrol.Earth Syst.Sci., 14, 765-781, 2010 www.hydrol-earth-syst-sci.net/14/765/2010/ (this corresponds to 1 kPa), and (3) the outflux plateau at the beginning of the experiment was not reproduced at all. Figure 15 shows a comparison between the accepted fit with the highest residuum and the rejected fit with the lowest residuum for all investigated soils.Using this procedure, the parameters were reasonably reproduced.The mean and standard deviation of the estimated values for each parameter is given in Table 4, the hydraulic functions for all converged fits together with the ones for the true parameters are shown in Fig. 16.The deviations in each column of Table 4 were calculated from the ensemble mean of converged results.When comparing the results of the inversion and its standard deviation calculated from the analysis of the sensitivity matrix with the real parameters and the variance encountered by the different fits in the ensemble, one can see that the standard deviations calculated from the analysis of the sensitivity matrix are too small in almost all cases.This is attributed to the fact that the sensitivity matrix does only give a linear approximation of the uncertainties.Deviations of the resulting fit parameters are quantified by the deviation coefficient d.No systematic deviations are found for the silt, however small systematic deviations are present for the sand and the sandy loam.This is expected, since the evaporation method is most sensitive in the dry range, where dynamics of the silt actually take place, but it is relatively insensitive for the wet range, the region of most of the dynamics of the sand and sandy loam soil.The silt has the smallest deviations since its dynamics are reaching most into the dry region of the three soils.This is also confirmed by the plots.For the water characteristic of the sand with 42 fits, 38 practically overlap while 4 deviate quite a bit.The sandy loam has the most significant deviations of all three soils in the water characteristic, only 12 of 27 fits are practically overlapping and 3 have small deviations, while the remaining 12 fits are relatively bad.The silt is nicely represented in all fits, only 2 of 34 fits are deviating slightly.Considering that the estimation of the hydraulic conductivity is difficult, the conductivity function is well represented for all soils.The analysis shows that applying the Monte-Carlo Levenberg-Marquardt approach for real experimental data would lead to a reasonable estimate of the parameters, since most fits converge to the correct parameters and the few deviating fits could be identified in the ensemble.

Conclusions
In accordance with experimental results from the literature, the physical analysis of the evaporation experiment model showed two different evaporation regimes.In regime I the outflux is limited by the diffusive soil-atmosphere boundary layer, while in regime II it is limited by the hydraulic properties.The sensitivity analysis showed that the sensitivity is especially high at the transition between the two regimes.In regime II, a drying front penetrates through the soil.At the location of this front, sensitivity is very high but this sensitivity peak can only partially be exploited as the potential drops below the measurement range of traditional tensiometers.However, if a profile measurement of water content is available the sensitivity peak allows to scan different layers of the sample.Positioning the tensiometer and the permittivity probe 2 cm below the sample surface gives a good sensitivity over the whole experiment.
As the transitions between these two regimes can be induced by boundary condition changes, a threestep experiment, where the water vapour pressure at the boundary is temporarily increased after regime II was reached, yields an increased sensitivity and a high measurement range.
The analysis of the response surfaces exhibited a single minimum for all parameters, which was mostly well localised.The estimation of the saturated hydraulic conductivity was improved significantly by adding a potential measurement to the outflux measurement.In contrast, adding a water content measurement yields no noticeable improvement.However, a combination of permittivity and tensiometer measurements gave nearly as good a response surface as a combination of evaporation flux and tensiometer measurements.However TDR probes which are typically used for water content measurements are no point measurements but have a rather high sampling volume, and the uncertainty for very dry conditions diverges.Therefore it is expected that for real measurements the response surface will get much worse.Furthermore uncertainties in porosity and soil matrix permittivity would introduce additional errors which were not considered in the simulations.Hence measuring j w will probably be much more accurate in real experiments.There is a significant improvement if the measurement range of the tensiometer is extended to −70 kPa.A hypothetical tensiometer with unlimited range shows a huge advantage.This emphasises the importance of extended-range tensiometers.
When excluding obviously diverged fits, the inverse model converged for the silt to correct solutions from all initial values.For the sand and sandy loam, some fits of the ensemble deviated significantly while others reasonably converged to the correct parameters.The deviations were attributed to the small sensitivity of the evaporation experiment in the wet range.Usage of the Monte-Carlo Levenberg-Marquardt approach allows a meaningful parameter estimation.Comparing the hydraulic functions estimated using the mean MCLM estimated parameters showed that only the water characteristic of the sandy loam deviated slightly and all others overlapped with the true functions.The distribution of the parameters caused by measurement errors and cross-correlation of the parameters was acceptable.

Fig. 1 .
Fig.1.Sketch of the experimental setup.Evaporation takes place into a gas-tight head space above the soil surface.Air is flowing through it to remove the water.Water vapour molar fraction and temperature in the head space is controlled to define the boundary condition.The water flux is measured by the difference in vapour content of the incoming and outgoing air in a controlled gas flow.

Fig. 2 .
Fig.2.Water flux at the upper boundary j w and potential in 2 cm depth of the onestep measurement.Notice that p w jumps from equilibrium to 1 kPa at t=0.

Fig. 6 .
Fig. 6.Sensitivity coefficients s(t) for the onestep experiment for the outflux (a), the potential in 2 cm (b) and the water content in 2 cm (c).A detail view of the first 100 h for the potential in 2 cm is shown in (d).Notice that p w jumps from equilibrium to 1 kPa at t=0.

Fig. 7 .
Fig. 7. Potential ψ m (z) (a), water content θ (z) (b) and hydraulic conductivity K(θ(z)) (c) distributions for the twostep experiment with 20 times higher α: (1) directly before the first boundary condition switch, (2) after relaxation, (3) after re-entering regime II and (4) at the end of the experiment.Notice the non-linear scaling of the depth axis.

Fig. 8 .
Fig.8.Influence of water vapour transport on the results: The threestep experiment was simulated taking vapour transport inside the soil into account (the normal case, red) and without considering it (blue).Notice that p w jumps from equilibrium to 0.25 kPa at t=0.For more explanations see Sect.3.2.1.

Fig. 9 .
Fig. 9.Sensitivity coefficients s(t) for the threestep experiment for the outflux (left), the potential in 2 cm (middle) and the water content in 2 cm (right).Notice that p w jumps from equilibrium to 0.25 kPa at t=0.

Fig. 10 .
Fig. 10.Relative sensitivity s(z,t i ) on α of the potential ψ m (top) and the water content θ (bottom), respectively, versus height z for different times t i , both for the threestep experiment.Note the switch in the sensitivity scale, marked with a dashed black line, and the non-linear scaling of the depth axis.

Fig
Fig. 11.Relative sensitivity s(z,t) of the potential ψ m on all variables at t=176.9 h (top).The potential and water content profor the undisturbed parameters as well as with α changed by 10 % (bottom) illustrate the source of the huge sensitivity peaks at z≈0.15 cm in the plot above, the location is marked with a dotted vertical line.Note the non-linear scaling of the depth axis.

Fig. 13 .
Fig.13.Same as Fig.12but with different target variable combinations.Note that real permittivity measurements have large errors in the dry range, which would lead to much worse results for the case c + ψ m .
. 12 but for different tensiometer tear-off values ψ limit with j w and ψ m as target variables.off values ψ limit with j w and ψ m as target variables.

Fig. 14 .
Fig.14.Same as Fig.12but for different tensiometer tear-off values ψ limit with j w and ψ m as target variables.
Fig. 15.Comparison between the accepted fits with the largest residuum and the rejected fits with the lowest residuum, for sand (a) (residuum 4900 (accepted) vs. 9700 (rejected)), loamy sand (b) (residuum 28 000 vs. 3300), and silt (c) (residuum 4600 vs. 7300), respectively.The reason for rejecting the curves were for (a) the deviation of the tensiometer, for (b) the first plateau was not represented, and for (c) the first peak was not represented, respectively.Note that for (b), although the potential of the accepted fit matches worse than in the accepted one, its deviation is just at the sort-out limit and limit cases were still accepted, while the plateau in the rejected fit is definitely not represented at all.

Table 2 .
The

Table 3 .
Uncertainties assumed for the virtual measurement devices.

Table 4 .
Results of the convergence study with j w and ψ m as target variables.For the resulting parameters and standard deviations from the model output, as well as the deviation coefficients d defined by Eq. (12), mean and variance of the inversion of the start parameter sets with each of the data sets are shown.