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.
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