Two of two topics about Evapotranspiration (and about the entropic origin of diffusion)

Did you care about the mass conservation, as expressed by the continuity equation ?

Evapotranspiration can be seen as a mass budget. Therefore, from this point of view, the above equation must hold. However we are not used to see it in its complete form. Usually in fact, the gradient of velocity is simply neglected. Or the other way around.  Whatever of the two, one could think that the right answer to the budget would come, obviously, to solve the Navier-Stokes (NS) equations for getting the velocity field to substitute inside the equation. 

But, this is not exactly the case. Precisely, to solve NS equations (numerically, clearly), we would have to put some boundary conditions at the bottom of the domain, i.e.  the water, the soil or the leaves, and this would be, hundred per cent, a no slip condition where velocity is null. As a matter of fact, if velocity would be null, we would have no flux of water vapour, implying that our deployment of the problem has something wrong. First consideration is that what we actually care is not the (mean) velocity of air but the mean velocity of water vapour itself and, while the (mean) velocity of air molecules at the boundary is null,  the mean velocity of water vapour molecules is not null cause to entropic forces, and set by the diffusion process. So, even at null (mean velocity), we have to introduce in our computing a Neumann (flux) condition based on some (non-equilibrium) thermodynamics.

If we are not going to solve appropriately the NS equation we can think to simplify the question substituting appropriate averages of the velocity field, for instance assuming that the flux is diffusive at the bottom (close to the surfaces) and logarithmic above, as derived from information of the Figure. This scheme has a long history and can be traced in the book by Brutsaert (1980). 

In either of the case, the interesting thing now is how to estimate the flux at the bottom from thermodynamics. A starting point for the investigations of these fluxes in the Einstein theory and the Teorell formula which express these fluxes as function of the chemical potential, and could be, eventually properly upscaled.

Still too vague .... I guess. But we can try to work on it.

One of two topics about Evapotranspiration (and the Penman Monteith equation)

Evapotranspiration (e.g. Brutsaert, Evaporation into the atmosphere,1980, or my lectures - in Italian - in English) is the phenomenon that describes collectively evaporation from water surfaces, and plants. Hydrologists used to used to neglect its estimation when they were concentrated on rainfall-runoff, and many, even nowadays, tends to forget its existence.

Once considered a stationary energy budget, as Penman (1948) showed, its estimation is reduced to a formula, the Penman formula, subsequently modified by Monteith (1965) to account for resistances due to soil suction and plants reaction to water stress (e.g. Rodriguez-Iturbe and Porporato, 2004).  The formula is universally knows as Penman-Monteith formula. Nobody exactly knows what the resistances  means, or better, how to estimate them, since they depends on many factors.  

However, being a formula, it does not require the solution of differential equations, and its evaluation depends only on the measure or the spatial extrapolation of the terms it contains. 

Radiation has actually its own difficulties to be estimated, but a reasonable way to do it, was found and improved during the years and is available (e.g. Formetta et al., 2013) .  Wind velocity can be either simulate or interpolate from measurements. On resistance terms either we can do educated guesses, or, after having estimated the potential evapotranspiration (PET), at least a first approximation can be found in appropriate factor of reduction linear proportional to the water content both for the soil cases (where the water content can vary between the residual water content and water content at saturation)  and vegetation case (where the extremes are, for each species, the wilting point and a critical water content, something below the air entry level). I would not have idea , I admit,  how to estimate  air humidity, if not interpolating measurements, but the air moist at saturation, thanks to the Clausius-Clapeyron relation depends only on temperature. 

For each of these quantities actually, in our ignorance we could give an estimate. If we further think to be able to assess the estimate on the measurements, we could think to be able to propagate di error to obtain (under the hypothesis of parameters distributed according to the Normal distribution) the "error" of our estimate. 

Starting from the formula in the figure, an from the hypothesis of normality, it is quite easy to find the distribution of evapotranspiration, which is actually itself a normal distribution whose mean is equal to the deterministic expression of ET, given, precisely, by formula in Figure above. 

So why don't try it ?

A more complex situation derives if we take into account space. In this case, any of the parameters in the game is a function of location, and the error of estimate too. However, using an appropriate linear model (for instance Kriging) is not overwhelming difficult to derive an estimate of both, and therefore coming out with a spatial estimate of ET and its error. 

Space clearly bring into the game the spatial variability of ET, evaporation rather than evaporation from soil or transpiration from plants (grass, shrubs and trees, or crops at least).  So, if the error of the parameters can be estimated, for estimating evapotranspiration itself further guesses need to be done about the spatial variability of soil cover and use. This, obviously,  brings inside other sources of uncertainty. However, at least some playing around with the appropriate tools (i.e. Jgrass-NewAGE) could be made. 

Obviously if stationarity of the energetic fluxes is not assumed, all it is another game (but you can use GEOtop, then).

A proposal of an R package for investigating some simple soil properties

In analysing the soil hydrological properties the first step is certainly to be able to draw the Soil Water Retention Curves (SWRC), at least according to the two major parameterisations: the van Genuchten and the Brooks and Corey ones. Possibly some outsider parameterisation could be also covered, as the double porosity type envisioned by Romano et al., and the (deprecated) Gardner ones (but used, for instance in TRIGRS to have analytical solution for the vertical infiltration). So why do not face this task to implement all of them in R (see also here, a previous post) ?
These parameterisations, in turn, could be used inside the Mualem scheme to obtain the hydraulic conductivity (K, here: slides no 127), and after plotting them (with ggplot2) their properties could be investigated. Particularly important it has been found to be the delay of increase/decrease with suction of K with respect SWRC (here, from slides no 72), so plotting them one above the other can be useful to understand qualitatively how the subsurface flow dynamics (excluding macropores) works at hillslope scale.  Plotting hydraulic capacity, C (slide 26), a.k.a. the derivative of the SWRC with respect to suction, could be nice to see. Also the plot of the hydraulic diffusivity,  D = K/C, can  be useful. So we have at least four functions for any parameterisation to plot.
One variation on this scheme can be considering the extension of the SWRC from the negative pressures (suction) domain to the positive ones. Since many phenomena occurs just at the edge of saturation a proper extension of the SWRC in this domain can help to understand. This paper by Schaap and van Genuchten can be used as a basis for any investigation.

To solve Richards equation analytically either adopt the Gardner parameterisation of SWRC or over-simplifiy it. Assuming constant hydraulic diffusivity (as a product of constant K and C) it is possible to obtain well know analytical solutions of the 1-D Richards equation derived from well-known heat transfer solution (here, from slides no 179 or here, in the original paper). These solution, not only offer an opportunity to plot another, more complicate function, but also allows for calculating an approximation of the 1-D infiltration phenomena induced by a variable rainfall, as a convolution of these solutions with the variable rainfall inputs (which requires a little more of R programming).
Once the whole machinery would implement, would be a piece of cake to evaluate the errors that can  be made with these solution by appropriately varying the parameters.

References

Baum R.L., Savage W.S., and Godt J.W., TRIGRS—A Fortran Program for Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Analysis, Version 2.0

D'Odorico, P., Fagherazzi, S., & Rigon, R. (2005). Potential for landsliding: Dependence on hyetograph characteristics. Journal of Geophysical Research, 110(F1), 1–10. doi:10.1029/2004JF000127

Rigon R. - Slides about Soil Hydrology (here in Italian, here in English), 2011

Rigon R. - An Overview of hillslope hydrology, 2013

Romano, N., Nasta, P., Severino, G., & Hopmans, J. W. (2011). Using Bimodal Lognormal Functions to Describe Soil Hydraulic Properties. Soil Science Society of America Journal, 75(2), 468. doi:10.2136/sssaj2010.0084

Schaap, M. G., and van Genuchten, M. T. (2006). A Modified Mualem–van Genuchten Formulation for Improved Description of the Hydraulic Conductivity Near Saturation. Vadose Zone Journal, 5, 27–34. doi:10.2136/vzj2005.0005

Richards 1D integration

As shown in detail at the July summer school about landslides, one of the first approximations to the full Richards equation is the 1-D  Richards equation (please take care of the fact that vertical is not slope-normal). If we want to solve it, besides using the celebrated Hydrus-1D, we can use the novel algorithm developed by Vincenzo Casulli and Paola Zanolli (finite elements, the first, finite differences the second - the second Authors actually implemented a 3-D algorithm). The conceptual line followed by the Casulli and Zanolli is "simple" (but simplicity, as usual, is the product of long and deep thinking):

•  discretise the system to obtain a non-linear system
•  use a Newton method to solve the non linear system

In turn, the Newton's method^1 is implemented in a way that it iterates over a sequence of linear equation in which the relevant equation matrix is symmetric, and therefore the Conjugate Gradient^2 (CG) method can be used to solve it efficiently.
That's the theory. And, in practice ?

S

You need to choose a language for implementing the algorithm. Assume you use Java, because you know it (or you have learnt it, or you just want to learn it, and this can be a challenging but nice exercise).  Since you will be doing a lot of matrix algebra, you first have to choose a library for solving linear algebra (or to implement it by yourself), And if the library can, without you very much take care of it,  thread your algorithms, it is a bargain.  One of these libraries is, for instance, the Parallel Colt. Obviously you need to get used to the library (after having installed it).  some other library can be more efficient, but in Parallel Colt you can find the Conjugate Gradient method already implemented. Once you manage both (PC and CG), you are pretty close to the solution of the problem (easy to say ;-) ) but ...

_________________________________________________________________________________
^1 - A good concise reference is the Davis' book "Direct methods for sparse linear system", SIAM 2006, but the appendix of the paper on Boussinesq equation can clarify several things.
^2 - Give a look to  An Introduction to the Conjugate Gradient Method Without the Agonizing Pain by Jonathan Richard Shewchuk.