Again on reproducible research

I've got this from a post on linked-in by Mariza Costa-Cabral, author of a pioneering paper on digital elevation models, and now senior scientists at Northwest Hydraulic Consultants, Inc.:

"I hope this is the start of something good for the sake of reproducibility of research: Mozilla analysis of scientific software. And if such automated or semi-automated debugging can be done, then I hope we can get to use it much BEFORE submitting the paper! :-) Imagine, how great it would be to get help with debugging! Also, please read the particularly sensible 2010 article by Nick Barnes: "

The paper had several comments around the web, including those by Titus Brown. Here further posts of mine on the same topic. 

SMAPEx and Rocco Panciera work in his seminar at Trento

Today Rocco gave his presentation here at his Alma Mater in Trento. Was a nice coming back and a beautiful presention of his work:

 Towards Global Remote Sensing of Soil Moisture: Australian Experiments

Abstract: Rocco panciera graduated from the University of trento in June 2013 with a thesis on distributed hydrological modelling under the supervision of Dr. Riccardo Rigon, and soon thereafter departed for the southern emisphere. He completed a Ph.D. at the University of Melbourne in 2010, with a thesis on remote sensing of near-surface soil moisture from airborne and spaceborne sensors, and since then has been working as a Research Fellow for the Australian Research Council, focusing on the estimation of soil moisture from airborne and spaceborne Synthetic Aperture Radar (SAR). Global monitoring of soil moisture using spaceborne passive microwave and SAR sensors is rapidly evolving from research into application, with the launch of soil moisture dedicated missions such as the Soil Moisture and ocean Salinity (SMOS, 2009) and the future Soil Moisture Active Passive (SMAP, 2014) missions. Moreover, there is an increasing trend in the availability of global coverage from Synthetic Aperture Radar (SAR) active microwave sensors such as SAOCOM, Sentinel-1, ALOS 2, Cosmo-SkyMed,  providing the opportunity for long-term temporally-dense series of microwave observations which are suitable for soil moisture monitoring at fine resolution (<1km).  Rocco was heavily involved in series of large-scale airborne experiments conducted in Australia in the 2005-2011 time frame that provided airborne data in support of algorithm development for soil moisture estimation from such missions. This presentation reports on such experiments and a number of research activities emerging from those related to the estimation of near-surface soil moisture and vegetation type and biomass from passive microwave and SAR observations.

Click on the figure to see the presentation slides.

Clarifications about a formula

I reading the TRIGRS manual I fond this:

Comparison of Computed Results with Published Analyses

We tested TRIGRS against Iverson’s (2000) examples from the Minor Creek landslide and from an experiment at the USGS debris-flow flume and were able to reproduce his results for pore pressure and factor of safety profiles. Thus, we were able to verify that the computer code accurately implements Iverson’s (2000) formulas for pore pressure and factor of safety. After publishing the original version of TRIGRS, students of Ricardo Rigon (University of Trento, Trento, Italy) pointed out an error in the coordinate transformation in Iverson’s (2000) original paper. Correction of that simple error causes our numerical results to differ somewhat from Iverson’s (2000) published values, though preserving their basic form. Folders named

25“MinorCreek” and “flume” contain files with necessary data to reproduce these results. By varying elapsed time, t , in the initialization file, the user can compute pressure head and factor of safety profiles like those of figures 7, 8, 10, and 11 of Iverson (2000). Note that the computed profiles conform to the limit imposed by equation 3. We have also verified that the new version of TRIGRS described in this report reproduces the published results of Srivastava and Yeh (1991). Data in the folder “sy91” produce profiles that allow the user to reproduce the curves plotted in figure 3 of Srivastava and Yeh (1991).

A part the misspelling of my name (two "r"s, not one like in Spanish), I need to precise that the error was discovered by Emanuele Cordano, a.k.a. the authors of at least  couple of papers with me, and some R packages useful for hydrologists.  His papers, maybe difficult to read, and so far not kissed by the goddess of citations,  remain among my own favorite. 

References

[ J28] - Cordano E. R. Rigon, A perturbative view on the subsurface water pressure response at hillslope scale, Water Resour. Res., Vol. 44, No. 5, W05407- W05407, doi:10.1029/2006WR005740, 2008

[J38] - Cordano, E., and R. Rigon (2013), A mass-conservative method for the integration of the two-dimensional groundwater (Boussinesq) equation, Water Resour. Res., 49, doi:10.1002/wrcr.20072.

Boussinesq code is available on Github.

R- Software

RMWAGEN, a weather generator, a package that contains functions for spatial multi-site stochastic generation of daily timeseries of temperature and precipitation. A presentation can be found here.

Soilwater, A package for plotting soil water retention curves and hydraulic conductivity written by Emanuele Cordano, Fabio Zottele and Daniele Andreis.

Luca Brocca seminar at Trento, September 25th, 2014

On september 25 Luca Brocca, one of the most promising young Italian hydrologists, will give a seminar on soil moisture modelling in  Hydrology. His presentation is already on-line on slideshare (click on the image below).

and therefore, you can have a preview of what he will be saying. Hope many will come. Living experience is, obviously, a must.
The abstract of the conference:

soil moisture governs the partitioning of mass and energy fluxes between the land surface and the atmosphere, thus playing a fundamental role for many scientific and operational applications, including flood forecasting, climate modelling, landslide prediction, numerical weather prediction, irrigation scheduling, to cite a few. Nowadays, soil moisture estimates from satellite sensors are becoming more readily available with a spatial and temporal resolution that is suitable for hydrological applications. Moreover, the accuracy of the satellite-derived soil moisture retrievals is found to be satisfactorily in many countries worldwide and mainly in the Mediterranean region.

This study aims at showing the reliability of the satellite soil moisture product derived by the Advanced Scatterometer (ASCAT) over Europe and, mainly, to understand how these observations impact the modelling of extremes in the Mediterranean region by considering three different applications.
The first application addresses flood forecasting (Brocca et al., 2010; 2012a), by assimilating the ASCAT soil moisture into a multi-layer continuous and distributed rainfall-runoff model, named MISDc. The Ensemble Kalman filter is adopted to optimally incorporate the soil moisture data into MISDc. Several catchments located in different climatic regions over Europe are used as case studies. Results reveal that the ASCAT soil moisture product can be conveniently used to improve runoff prediction, mainly if the soil wetness conditions before a storm event are highly uncertain or unknown. However, reliability differs according to the climatic region, the soil/land use conditions and the size of the catchments under investigation. Therefore, the open issues that should be addressed in future studies are also given.

The second application investigates the use of the ASCAT soil moisture product for predicting the movement of a rock slope located in central Italy, the Torgiovannetto landslide (Brocca et al., 2012b). By using a statistical approach, the opening of the tension cracks, recorded by an extensometers network operating in the area, as a function of rainfall and soil moisture conditions prior the occurrence of rainfall, are predicted in the period 2007-2009. Results indicate that the regression performance (in terms of correlation coefficient) significantly increases if the ASCAT soil moisture product is included.

Finally, the third application aims at estimating rainfall starting from soil moisture observations (Brocca et al., 2012c). Specifically, by inverting the soil water balance equation, a simple analytical relationship for estimating rainfall accumulations from the knowledge of soil moisture time series is obtained. Satellite and soil moisture observations from three sites in Europe are used to test the developed approach that showed reasonable results thus opening new opportunities for rainfall estimation at catchment/global scale. 

OCNs are back!

OCNs stands for optimal channel networks. They are topological configurations of channels that minimize energy expenditure and reproduce the statistical features of real river networks. They were the main topic of my Ph. Thesis, and go some influence in understanding why in nature, many structures are networks, and two of my most cited papers cover this topic.  Eventually my advisors and co-workers pursue further the studies and pushed forward the research.

In figures above (top) a small random network developed according to an Eden growth procedure which start from the outlets (seeds of the growth), and (bottom) a OCN developed from it.
Recently I was asked to  go back to the topic and do some simulation. For the task, I decided to abandon the old version of the program (first programmed in Mathematica, and secondly in C) and to build a new one in Java. The new version is part of the Java for Hydrologists (and geoscientists) 101 and can be found at github under the project Geoframe (package OCNs). News will follow.

Worth to say, other groups of people are working on OCNs. Notably Lily Briggs and her advisor, prof. Mukkai Krishnamoorthy, are also working on OCNs. She built three-dimensional OCNs. Her paper is  here and her code is here, in Github.

Earth System Modelling Framework and Related Resources

Well, we pushed on OMS and Java. However, for FORTRANists (or FORTRANers ?) a very interesting and modern resources is the Earth Science Modelling Framework or (ESMF in short). This was due to a large effort by NOAA, UCAR and others. Learning it has an apparent steep learning curve.  In any case, either you use it or do not use it, it is certainly a term of reference. Especially for the tools and the concepts developed, independently from their actual implementation. Click on the image to get the sites.

I obviously knew the program since a lot of time. However, It was brought to my attention again in these days because I was looking for the models' run metadata implementation, which is now part of the Earth System Documentation effort. The general goal of the initiative can be found in this 2012 Fall AGU Meeting  presentation.

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

The JGrass-NewAGE informatics

As I tried to convey in previous posts, since more than five years, I am working to the idea to built a hydrological model by components. Well, this should have been the first paper in row, but as often happens, is the fourth to have been submitted.

The paper contains, by Formetta et al., the main ideas behind this type of modelling and shows that the system actually works, it is not just a matter of  speculations. We did it. At the moment it is at a late stage of review on Environmental Modelling & Software, and we eventually ask for it being open access.

The rational of the paper is expressed in its introduction:

"Scientists demand more and more the availability of simulation model’s source code since it has become a key factor for the understanding, validation, and advancement of science (e.g. Ince et al. [2012]). However, this is not enough, even if the source code would be available, the growing complexity of modelling code makes model development progress challenging to understand and manage. In fact, if model code distribution is matter of policies (e.g. Annan et al., [2013]), external inspection and analysis of models, improvement and contribution are difficult or even impossible when the software is inadequately engineered. The implementation of many environmental processes intimately interlinked (as snow, runoff production, evapotranspiration, in the hydrological case), is usually difficult to understand per se, but models writing in traditional monolithic forms, as defined in Rizzoli et al. [2005], makes their implementations overwhelming hard to follow, and models themselves practically impossible to be verified (Quesnel et al. [2009]). As a matter of facts, the traditional modelling practices preclude easy understanding, rapid reuse and improvement of the source code, and eventually obstacle seamless advancement in science.

Part of models’ obscurity has its foundation in bad documenting practices, and many researchers’ community that rely on computational methods and techniques as part of their day-to-day activities proposed shared protocols, like the so called Overview, Data concepts, Details method (ODD, e.g. Grimm et al., [2006]), to improve documentation effectiveness. However, “reproducible-research systems” (RRSs), making easier to document any step during research from model implementation and data preparation to output analyses, would greatly help correct policies to be adopted. Actually from a RRS system we would expect, besides model codes sharing, tools to allow the researchers, on one hand, to repeat the simulations in the same conditions and, on the other hand, to spend more time on their own science.

Many of software infrastructures or modeling frameworks (MF) were actually designed and built to streamline the process of a sound scientific production (e.g. Wesselung et al. [1996], Argent [2005], Rizzoli et al. [2005]). Among those that specifically target the support of hydrological modelling are the Spatial Modelling Environment (SME, Maxwell and Costanza [1997]), The Invisible Modelling Environment (TIME, http://www.toolkit.net.au/Tools/TIME) and hydrological derivative tools like, E2 (Argent et al. [1999]), OpenMI (http://www.openmi.org/), Moore and Tindall [2005]), and the Object Modelling System (OMS, David et al. [2002, 2013]), Common Component Architecture (CCA, Bramley et al. [2000]) and Earth System Modeling Framework (ESMF, Hill et al. [2004]).

However, most of the above MF require a quite significant learning curve that not all scientists, even proficient modellers, are willing to make.
Therefore, in order to ease the transition into modern programming environments, some modelling efforts and projects recently focused on providing code generation support and reducing the invasiveness of frameworks (Lloyd, 2010) into the model. Especially the third version of OMS and the BioMA project (BioMA, 2012) revealing promising perspectives.

A RSS would not be complete without including data visualization. Gardner and Manduchi (2007), among others, emphasize that in order to optimize scientific productivity, a RRS infrastructure should include not only the computational cores but also visualization and data-processing tools necessary to synthesize knowledge from high volumes of inputs and outputs.

Indeed, tools of choice for the visualization of hydrological processes have been for a long time Geographic Information Systems (GIS) (Maidment [1993]; Grayson et al. [1992]). However, traditional GIS are usually designed for managing static, non- temporal information layers. They are not designed to interact with the dynamic modelling (e.g. Burrough et al. [1998]; Wesselung et al. [1996]). In fact, the interaction between models and GIS can be described as “off-line” and it is performed with integration strategies that affect either the functionality of GIS tools or the usability of models.

Instead, the MF listed above offer a proper abstraction to streamline the interaction with a GIS. They promote the separation of the model into well defined modules or components, each with a well-defined way to interact with others through specified interfaces. Through their interfaces the modules can communicate and exchange data.

Therefore it is also timely for GIS and hydrological model components to constitute a pool of interoperable tools that can be blended together to create software that is accurately tailored to geosciences."

The revised version of the paper prior to publication, in pdf format, can be retrieved from here (or clicking below the picture). The final version is here instead.