Key physical factors and processes affecting non-uniform and preferential flows in porous media

Preferential and non-uniform flows are induced by biotic (e.g., earthworms and roots) and abiotic factors and processes (e.g., wet-dry and freeze-thaw cycles, lithology, and structure) as well as anthropogenic activities (e.g., tillage and cultivation methods in agricultural land, human-made landforms from waste rock dumping and disposal strategies in mining). The understanding of preferential flow is of premium importance in relation to soil hydrology, as it can move a considerable amount of water and solutes in porous media. Preferential flows can occur spatially from the pore scale to entire catchments, across large regions. Temporally, the preferential processes can change during hydrological events, from within hours to seasonal events, and across inter-annual variations of years.
This session welcomes studies on experimental and theoretical challenges to identify, quantify, and model the key physical factors and processes that are responsible for preferential flows in porous media across scales (from pore scale to catchment scale). Contributions are welcome to reflect on experimental studies, novel approaches, and advances in solutions to:
• Understand the geometry and connectivity, formation and dynamics of fissures, fractures, and macropores and their effect on preferential flow;
• Understand the effect of physical processes and geochemical processes on the dynamics of macropores and fracture networks;
• Develop and refine models for quantifying preferential flow, from pore scale to pedon scale and entire catchments and landscapes;
• Unpacking the pore structure of soil using new methods and approaches, including the use of non-Newtonian fluids, for improved characterization of heterogeneous soils and the advancement of flow and transport modeling.

Convener: Laurent Lassabatere | Co-conveners: Majdi R. Abou Najm, Rafael Angulo-Jaramillo, Thomas Baumgartl, Jannes KordillaECSECS, Mandana ShayganECSECS
vPICO presentations
| Wed, 28 Apr, 13:30–15:00 (CEST)

vPICO presentations: Wed, 28 Apr

Solicited talk
Majdi R. Abou Najm and Keith Beven

Peter Germann died on December 6th 2020 in Bern, Switzerland. Known for a wide range of contributions to the physics of soil-water interactions and flow, his name (along with Keith Beven, his career-long collaborator and fiend) is recognized by an entire generation of soil physicists and hydrologists who studies macropore and preferential flows. They both co-authored the classic, and highly cited 1982 review paper in Water Resources Research on Macropores and Water Flow in Soils. Peter’s PhD work between 1976-1980 was a study of soil-water relations based on maintaining a network of 35 nests of tensiometers at 10 different depths down to 3m. At that time, these were still manual tensiometers coupled to mercury manometers that were read every 2 to 3 days for 3 years. One of the features that this remarkable data set revealed was that during infiltration, wetting in some cases occurred at depths, apparently by-passing the tensiometers above. This is what we all now know as preferential flow. Another was the large heterogeneity in responses between sites and between wetting events. For the major part of his research career, Peter was a strong advocate for a reconsideration of the physics of water flow through soils and, in particular, for the limitations of the Darcy-Buckingham-Richards flow theory. Peter later developed the kinematic wave approach into a theory of viscosity (rather than capillarity) dominated film flows subject to Stokes’ law during infiltration. He summarised his research work in his 2013 book on the subject published by the University of Bern. Peter held academic positions at the University of Virginia in Charlottesville, at Rutgers University, and at the University of Bern back in Switzerland where he stayed until he retired in 2009, and held an Emeritus position until 2015.   He continued to publish papers until shortly before his death which followed 2 major strokes. In this talk, we will go over Peter’s main contribution and research highlights in the area of macropores and preferential flows. Peter was no stranger to EGU, and many know him and have met him in this session or others. For those who knew Peter, they will miss his enthusiasm, his critical mind, his genuine care for the state of soil physics, his thoughtful responses, and his humour. He was a great source of inspiration to us and many others. Peter will be missed by many in soil science.

How to cite: Abou Najm, M. R. and Beven, K.: Preferential flows and the legacy of Peter F. Germann (1944-2020), EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9412,, 2021.

Preferential flow and infiltration experiments
Felix Abayomi Ogunmokun and Rony Wallach

Preferential flow pathways and uneven soil water and chemical distribution are intrinsic phenomena in water repellent soils. These uneven water and chemical distribution reduce water uptake by the plant roots on one hand and enhance deep percolation and chemical leaching, on the other hand, thereby enhancing soil and groundwater pollution. The results of attempts to remediate soil water repellency and heterogeneous spatial distribution of soil moisture and chemicals within the root zone by surfactant application will be addressed.  

This study was conducted in a commercial citrus orchard in central Israel that is irrigated with treated wastewater. Previous studies have revealed that prolonged irrigation using treated wastewater renders the soil water repellent with its associated adverse effects. The soil water distribution within the soil profile was monitored by frequent electrical resistance tomography (ERT) scans. The spatial distribution of different chemicals within the soil profile was obtained by chemical analysis of disturbed soil samples taken manually along a line transects. Two methods of surfactant application were used and compared: 1) on soil surface spraying (area source), 2) via drippers application (point source).

Surfactant spraying onto the water repellent soil's surface succeeded in turning the soil wettable, diminishing the preferential flow pathways, and renders the soil water and dissolved chemicals uniformly distributed. In contrast, drip applied surfactant exacerbated the incidence of preferential flow pathways and the leaching of solutes from the soil. Moreover, the overall average water content in the 0-40 cm soil layer significantly increased with surfactant spraying than with drip application even though both were higher than the control plots. These results substantiate previous laboratory-scale studies in which surfactant was applied to water repellent soils packed in a transparent flow chamber by these two methods. Additionally, the yield from the on-surface surfactant sprayed plots show a slight continuous increase compared to the untreated plots.

How to cite: Ogunmokun, F. A. and Wallach, R.: Remediation of intrinsic preferential flow pathways in water repellent soils' profile using surfactant, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8537,, 2021.

Ryan Stewart, Majdi R. Abou Najm, Simone Di Prima, and Laurent Lassabatere

Water repellency occurs in soils under a wide spectrum of conditions. Soil water repellency can originate from the deposition of resinous materials and exudates from vegetation, vaporization and condensation of organic compounds during fires, or the presence of anthropogenic-derived chemicals like petroleum products, wastewater or other urban contaminants. Its effects on soils range from mild to severe, and it often leads to hydrophobic conditions that can significantly impact the infiltration response with effects extending to the watershed-scale. Those effects are often time-dependent, making it a challenge to simulate infiltration behaviors of water-repellent soils using standard infiltration models. Here, we introduce a single rate-constant parameter (αWR) and propose a simple correction term (1-e-αWRt) to modify models for infiltration rate. This term starts with a value of zero at the beginning of the infiltration experiment (t = 0) and asymptotically approaches 1 as time increases, thus simulating a decreasing effect of soil water repellency through time. The correction term can be added to any infiltration model (one- two- or three-dimensional) and will account for the water repellency effect. Results from 165 infiltration experiments from different ecosystems and wide range of water repellency effects validated the effectiveness of this simple method to characterize water repellency in infiltration models. Tested with the simple two-term infiltration equation developed by Philip, we obtained consistent and substantial error reductions, particularly for more repellent soils. Furthermore, results revealed that soils that were burned during a wildfire had smaller αWR values compared to unburned controls, thus indicating that the magnitude of αWR may have a physical basis.

How to cite: Stewart, R., Abou Najm, M. R., Di Prima, S., and Lassabatere, L.: A quick fix for modeling infiltration in water-repellent soils, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-3408,, 2021.

David Moret-Fernández, Borja Latorre, Laurent Lassabatere, Simone Di Prima, Mirko Castellni, Deniz Yilmaz, and Rafael Angulo-Jaramillo

The 3-D Haverkamp et al. (1994) model for disc infiltrometer measures on homogeneous media involves the following parameters: the soil sorptivity, S, the saturated hydraulic conductivity, Ks, the β parameter and the A= (γ S2)/(rd*Δθ) term, where rd is the disc radius, Δθ is the soil water increase and γ is proportionality constant. Fixed β and A values are commonly used in most cases. S, and Ks can be estimated from the inverse analysis of a cumulative infiltration curve by fitting it the Haverkamp model. For practical reasons, Haverkamp implicit model is replaced by its 4-term (4T) approximate expansion for the transient state. The first part of this work analyzes the influence of layered soils on Ks and S estimates, and designs a new procedure, sequential Analysis of Infiltration curve (SAI), for treating infiltration curves impacted by soil layering. The SAI method analyzes a sequence of increasing dataset for a given infiltration curve and fits to the 4T expansions to estimate Ks, S. Then estimates and RMSE are reported as a function of the number of data points used for the fit. The method was applied on synthetic profiles with homogeneous loam soil, six layered profiles involving a 1, 2 and 3 cm thickness loam layer over silty or sandy loam soils, respectively. Erroneous estimates of Ks and S were obtained when the total infiltration curves were considered for the analysis, regardless of the presence of soil layering. In opposite, estimates were improved using the SAI method for the layered systems. The SAI method relies on the fact that the RMSE increases when the wetting front reaches the interface between the upper layer and the lower layer. Such increase allows (i) the detection of the soil heterogeneity, (ii) the determination of the optimum infiltration time, to, that corresponds to the minimum value of RMSE, and, (iii) accurate estimation the upper layer Ks and S.

Taking use of the SIA procedure, the second part of this communication studied the relationship between β and A, and proposed a new procedure to improve the estimate of Ks and S and approach β. The analysis was applied on synthetic infiltration curves simulated on homogneneous and layered columns. The results showed that different combinations of β and A resulted in similar Ks. Overall, optimization of Ks, S and A for different β values showed that β had an important effect on A and Ks, but not on S and RMSE.  We propose approaching the optimum β as the β for which is closer to zero, where A and Aexp are the optimized and measurable parameter, respectively. While the optimum β is calculated, Ks and S are computed by applying the optimum β to the respective quadratic β(Ks) and β(S) relationships. This methodology allowed improving the estimate of Ks giving good approaches of β (36% error) and omitting the erroneous praxis of using constant β and A values.

Haverkamp, R., et al. 1994. 3. Water Resources Research 30, 2931–2935.

How to cite: Moret-Fernández, D., Latorre, B., Lassabatere, L., Di Prima, S., Castellni, M., Yilmaz, D., and Angulo-Jaramillo, R.: Analysis of 3D infiltration curves measured with disc infiltrometer in heterogeneous soil profiles: Sequential analysis of infiltration data and estimate of β, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12442,, 2021.

Deniz Yılmaz, Laurent Lassabatere, David Moret-Fernández, and Borja Latorre

For heterogeneous soils, accurate water modeling in unsaturated soil conditions is a very important prerequisite since activation of macropore during the flow process is directly linked to the bulk saturation of the soil matrix. Indeed, macropores activate and begin to infiltrate when they receive the runoff from saturated matrices. To this point, the accurate estimation of the matrix hydraulic properties is of uttermost importance. We then focus on the accuracy of estimates for hydraulic parameters, by fitting to the well-know Haverkamp 1D analytical infiltration equation, that is widely used for Beerkan type infiltration. This equation involves an infiltration constant called β that is fixed to a by-default value of 0.6. This value is considered for relatively dry condition and for all type of soils, including fine matrices (silt, clay, etc.) but also coarse soils which are prone to preferential flows. However, the values of β have already been questioned by several authors. In this study, we performed a numerical study to investigate the value of β. Several cumulative infiltrations were numerically generated and fitted to Haverkamp’s model to derive the parameter β. This was then plotted as a function of initial water content and the type of soil. We proved that β is not constant. Especially for lower permeable soils, previous studies point that β value must be over 1 which is in contradiction with the domain of definition of β and the usual ranges considered for this parameter. Therefore, using β equal to 0.6 leads to an overestimation of Ks, leading to an overestimation of the soil capability to infiltrate and the prediction of the water budget. Numerical investigations of β show that this parameter is also a function of the degree of saturation. As defined by Haverkamp (1994), it varies from 1 for dry soil conditions to zero for saturated conditions. The hypothesis of the constancy of β allows easy integration of Richards 1D equation, leading to the formulation proposed by Haverkamp et al. However, it conducts for low permeable of soils to overestimate Ks. In this study, we demonstrate that the use of the adequate function for describing β in function to the degree saturation and the soil type improves significantly the accuracy of Haverkamp’s model.

How to cite: Yılmaz, D., Lassabatere, L., Moret-Fernández, D., and Latorre, B.: Investigation of β parameter used in Haverkamp 1D-analytical infiltration equation , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12851,, 2021.

Laurent Lassabatere, Simone Di Prima, Paola Concialdi, Majdi Abou Najm, Ryan D. Stewart, Vincenzo Bagarello, Massimo Iovino, Mirko Castellini, Jesús Fernández-Gálvez, Joseph Pollacco, Deniz Yilmaz, and Rafael Angulo-Jaramillo

Preferential flow is more the rule than the exception. Water infiltration is often led by preferential flow due to macropores, specific soil structures (e.g., aggregates, macropore networks), or lithological heterogeneity (occurrence of materials with contrasting hydraulic properties). Water infiltration in soils prone to preferential flow strongly depends on soil features below the soil surface, but also the initiation of water infiltration at the surface. When the macropore networks are not dense, with only a few macropores intercepting the soil surface, water infiltration experiments with ring size in the order of 10-15 cm diameter may overlook sampling macropore networks during some infiltration runs, minimizing the effect of macropore flow on the bulk water infiltration at the plot scale.

In this study, we investigated the effect of ring size on water infiltration into soils prone to preferential flow. We used two ring sizes: small (15 cm in diameter) and large (50 cm in diameter). By doing so, we hypothesized that the large rings, sampling a more representative soil volume, will maximize the probability to intercept and activate a macropore network. In contrast, the small rings may activate the macropore network only occasionally, with other infiltration runs mainly sampling the soil matrix. Thus, the small rings are expected to provide more variable results. On the other hand, the large rings are expected to provide more homogeneous results in line with the soil's bulk infiltration capability, including all pore networks at the plot scale.

Three different sites were sampled with varying types of preferential flow (macropore-induced versus lithological heterogeneity induced). The experimental plan included inserting large rings at several places in the experimental sites with a dozen small rings nearby to sample the same soil. All the rings were submitted to a similar positive constant water pressure head at the soil surface. The cumulative infiltrations were then monitored and treated with BEST algorithms to get the efficient hydraulic parameters. Note that the cumulative infiltration could not be compared directly since lateral water fluxes varied in extent and geometry between the different ring sizes. The impacts of the ring size on the magnitude of cumulative infiltration and related estimated hydraulic parameters were discussed. Our results demonstrated the impact of ring size but also the dependency of such effect on the site and the type of flow.

Our results contribute to understanding preferential flow in heterogeneous soils and the complexity of its measure using regular water infiltration devices and protocols.

How to cite: Lassabatere, L., Di Prima, S., Concialdi, P., Abou Najm, M., Stewart, R. D., Bagarello, V., Iovino, M., Castellini, M., Fernández-Gálvez, J., Pollacco, J., Yilmaz, D., and Angulo-Jaramillo, R.: Coupling large and small ring infiltration experiments for investigating preferential flow, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-7355,, 2021.

Effects of preferential flow on reactive transfer
Ekaterina Kolchanova and Nikolay Kolchanov

We study convective instability in the vertically layered porous media saturated with mixture. The mixture consists of a carrier fluid and solid nanoparticles. The nanoparticles are considered as solute within the continuous approach. The porous media are two horizontal sublayers with different permeabilities. The solute concentration is maximal near the upper boundary and is zero near the lower boundary of the superposed sublayers. Thus, one has suitable conditions for the onset of solutal convection in the gravitational field.

The porous sublayers are reactive media, which can absorb nanoparticles. The mixture transport here is accompanied by immobilization. It is described by the mobile/immobile media model. The mobile phase is carried by fluid flow, while the immobile phase is absorbed by porous matrix. The linear kinetic equation for the mixture redistribution between the phases is applied. The Boussinesq approximation is used in the equations for convection in each of the sublayers. Numerical simulation is performed by the shooting method.

We apply a linear stability theory to find the threshold Rayleigh-Darcy number for the onset of solutal convection. This similarity criterion is determined through the average permeability and porosity of uncontaminated porous sublayers. For the first time, we introduce a solutal pore shrinkage coefficient, which is analogous to the thermal expansion coefficient for thermal natural convection. This coefficient shows that porosity decreases as the concentration of immobile phase grows in the presence of sorption. Particles in this case join the surface of pores and shrink the void space.

Firstly, we find the permeability ratios for bimodal marginal stability curves in the uncontaminated sublayers. Here, the sublayer permeabilities differ by approximately 100 times. The bimodal curves demonstrate the competition between two convective instabilities. One of them is for the local convective rolls that generate within the more permeable layer and the other is for the large-scale rolls penetrating both layers. The rolls are similar to thermal natural convection in the multi-layered porous media studied by McKibbin and O'Sullivan (1980). For sorbing porous media, the type of convective rolls strongly depends on the solutal pore shrinkage coefficient. Even a small change in its value can produce a large variation of flow streamlines from the convective rolls localized within the upper highly permeable sublayer to the rolls covering both the upper and lower sublayers. The observed sorption effect on the transition from local to large-scale convection is due to the fact that the permeability ratio depends on the solutal pore shrinkage coefficient. It is also found that sorption effect delays the onset of solutal convection.

The work was supported by the Russian Science Foundation (Grant No. 20-11-20125).

How to cite: Kolchanova, E. and Kolchanov, N.: Solutal Convection in Layered Sorbing Porous Media, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-3852,, 2021.

Luis Alfredo Pires Barbosa and Horst H. Gerke

During preferential flow events, soil macropores such as cracks and biopores (decayed root channels and earthworm burrows) may allow water and solutes to bypass the lower permeable soil matrix. The biopore walls are inherently compacted by the locomotion mechanism employed by earthworms and roots. In addition, there are the excretion of biopolymers and the hydrophobicity of mucilage excreted by the roots of plants or mucus by earthworms. This gives to the biopore a coating with physicochemical properties distinct from the soil matrix, such as wettability, sorption and cation exchange capacity. Consequently, changes in the mechanical properties of the material in that region are also expected, ensuring greater mechanical stability.

However, micro structural features (i.e. crack size distribution) are still poorly explored. The objective is to analyze such features in detail, in order to better understand the effects of the coating material on soil macro mechanical behaviour (i.e. tensile strength) to explain the flow exchange between biopore and the soil matrix.

Therefore, soil samples were collected from Bt horizons of two Haplic Luvisols located in northern Bohemia (Hnevceves, near Hradec Kralove, Czech Republic; 50°18′47′′ N, 15°43′03′′ E). From these air dried samples, three earthworm burrows were identified and carefully separated from soil matrix.

The samples were scanned with X-ray microtomography (X-TEk XCT 225, Nikon Metrology), using 100 keV, 120μA and no filter. The reconstruction of three-dimensional images was done with the CT Pro 3D software package (version 3.1) at a spatial resolution of 10μm and 8-bit gray scale resolution. The permeability in each region was calculated along the biopore and perpendicularly to the biopore from matrix to coating using stokes solver.

The calculated hydraulic permeability for coating and matrix was 55 and 0.4 μm2 along biopore direction and 11 and 3 μm2 perpendicularly to the biopore. The results from image analysis show no differences in crack size distribution between the materials, but the number of cracks and connections were superior for the coating material, suggesting that the differences in the pore structure can strongly affect the macropore-matrix mass exchange.

How to cite: Pires Barbosa, L. A. and Gerke, H. H.: Structural characterization and permeability estimation for soil regions at the interface between biopore coating and soil matrix, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15959,, 2021.

Lian Zhou, Laurent Lassabatere, Jean-François Boily, and Khalil Hanna

Mass transport is significantly impacted by the nature of flow and, in particular, the occurrence of preferential flows. Most of the time, studies focus on observing preferential flow and its impact on mass transport either at the lab or the field scales. In the lab, real matrices are considered and embedded into columns, and mass transport is assessed for specific solutes and under controlled conditions (constant flow rate, saturation degree, etc…). However, very few studies use synthetic matrices and need to face matrix complexity in terms of both physics and chemistry. Such a complexity provides noise, uncertainty, and difficulty for the clear identification of mechanisms. This study made use of synthetized goethite nanoparticles as the reactant (sorption sites) combined to standardized sand to make a synthetic well-controlled porous medium. The goethite texture was changed during its fabrication to form two types of goethite-sand mixture: goethite-coated sand and goethite-aggregated sand. In the first case, goethite particles deposit at the surface of sand grains (forming a kind of coating), whereas goethite forms aggregates in the second case. The two types of columns were submitted to the injection of a tracer and two solutes: nalidixic acid (NA) and silicate. Our results show that flow remains mostly homogeneous, with the tracer following a straightforward ADE advection Dispersion Equation) process and no water fractionation into mobile and immobile water fractions. The minimal content of goethite (in the order of a few percent) does not change flow pathways. In contrast, the reactive transfer of NA and silicate is significantly impacted with less sorption, and much more solute spread in goethite-aggregated columns. NA and silicate cannot reach sites inside aggregates, reducing and slowing down their adsorption. In other words, changing the deposition mode of goethite nanoparticles on sand did not impact most of the flow and non-reactive transfer. It however greatly impacted reactive transfer. In addition, our results show that even if tracer experiments are performed for columns and attest of homogeneous flow, great care must be taken for reactive solutes. Tracers may not be the right tool to provide a clear picture of local hydraulic conditions at the vicinity of sorption sites, which are of utter importance for understanding reactive solute transfer.

How to cite: Zhou, L., Lassabatere, L., Boily, J.-F., and Hanna, K.: Conditioning preferential flow by using synthetic goethite aggregates, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-7429,, 2021.

Jérôme Raimbault, Laurent Lassabatere, Pierre-Emmanuel Peyneau, Denis Courtier-Murias, and Béatrice Béchet

Preferential flow is quite usual in natural environments. Non-uniform and preferential flows co-exist or alternate, impacting water transport and contaminant transfer through the vadose zone. In this study, we investigated how macropore-induced flow affects manufactured nanoparticles, as emerging contaminants reactive transfer. Previous studies showed that the presence of a macropore into water-saturated soil columns can foster preferential water flow within the macropore. One could expect that this preferential flow may increase contaminant transfer and reduce retention by the matrix in the case of contaminant, as previously reported. In this study, we injected pulses of silver nanoparticles to assess their transfer through sand columns with and without a macropore. Both systems (with and without macropore) were studied under similar conditions. An unexpected result was obtained: more nanoparticles were retained in the system with a macropore, i.e., with a preferential flow. This result is quite counter-intuitive. It appears that the relation between flow homogeneity and contaminant retention is not straightforward. Some possible explanations, related to chemical and physical kinetics, are put forward to explain the experimental results.

How to cite: Raimbault, J., Lassabatere, L., Peyneau, P.-E., Courtier-Murias, D., and Béchet, B.: May macropores increase contaminant retention? , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-10371,, 2021.

Modelling and integrative approaches
Swamini Khurana, Falk Heße, and Martin Thullner

In a changing climate scenario, we expect weather event patterns to change, both in frequency and in intensity. The subsequent impacts of these changing patterns on ecosystem functions are of great interest. Water quality particularly is critical due to public health concerns. Already, seasonal variation of water quality has been attributed to varying microbial community assemblages and nutrient loading in the corresponding water body but the contribution of the variations in the quantity of groundwater recharge is a missing link. It is thus beneficial to establish links between external forcing such as changing infiltration rate or recharge on nutrient cycling in the subsurface. We undertake this study to investigate the impact of temporal variation in external forcing on the biogeochemical potential of spatially heterogeneous subsurface systems using a numerical modeling approach. We used geostatistical tools to generate spatial random fields by considering difference combinations of the variance in the log conductivity field and the anisotropy of the domain. Tuning these two parameters assists in effective representation of a wide variety of geologic materials with varying intensity of preferential flow paths in the heterogeneous domain. We ran simulations using OGS#BRNS that enables us to combine a flexibly defined microbial mediated reaction network with the mentioned spatially heterogeneous domains in transient conditions. We propose that a combination of estimated field indicators of Damköhler number, Peclet number (transformed Damköhler number: Dat), and projected temporal dynamics in surface conditions can assist us in predicting the change in biogeochemical potential of the subsurface system. Preliminary results indicate that we miss potentially critical variations in reactive species concentration if we neglect spatio-temporal heterogeneities for regimes where 1<Dat<40. For regimes characterized by values outside this range, we propose that spatio-temporal heterogeneities due to subsurface structure and changing hydrological forcing may not be relevant.

How to cite: Khurana, S., Heße, F., and Thullner, M.: Predicting change in biogeochemical potential of subsurface systems with changing hydrogeological conditions, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-6349,, 2021.

Robert Mietrach, Thomas Wöhling, and Niels Schütze

The classical formulation of Richards' equation is relying on a unique functional relationship between water content, conductivity and pressure head. Some phenomena like hystersis effects in the water content during wetting and drying cycles and hydraulic non-equillibrium cannot be accounted for with this formulation. Therefor it has been extended in different ways in the past to be able to include these effects in the simulation. Each modification comes with its own challenges regarding implementation and numerical stability.
The Method Of Lines approach to solving the Richards' equation has already be shown to be an efficient and stable alternative to established solution methods, such as low-order finite difference and finite element methods applied to the mixed form of Richards' equation.
In this work a slightly modified Method Of Lines approach is used to solve the pressure based 1D Richards' equation. A finite differencing scheme is applied to the spatial derivative and the resulting system of ordinary differential equations is reformulated as differential-algebraic system of equations. The open-source code IDAS from the Sundials suite is used to solve the DAE system. Different extensions to Richards' equation have been incorporated into the model to address the shortcomings mentioned above. These extensions are a model able to simulate preferential flow using a coupled two domain approach, a simple hysteretic model to account for hysteresis in the water retention curve and also two models to either fully or partially calculate hydraulic non-equillibrium effects. To verify the numerical robustness of the extended model, stochastic parameterizations were generated that represent the full range of all soil types. Simulations were carried out using these parameter sets and real-world meteorological boundary conditions at 10 minutes time intervals, that exhibit drastic flux changes and poses numerical challenges for classical solution methods.

The results show that not only does the extended model converge for all parameterizations, but that numerical robustness and performance is maintained. Where applicable the results have been verified against solutions from the software Hydrus and show good agreement with those.

How to cite: Mietrach, R., Wöhling, T., and Schütze, N.: Modelling non-equillibrium unsaturated flow in soils during sudden pressure head changes by solving Richards' equation with a Method of lines approach, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15546,, 2021.

Jesús Fernández-Gálvez, Joseph Pollacco, Stephen McNeill, Sam Carrick, Linda Lilburne, Laurent Lassabatere, and Rafael Angulo-Jaramillo

Hydrological models use soil hydraulic parameters to describe the storage and transmission of water in soils. Hydraulic parameters define the water retention, θ(ψ), and the hydraulic conductivity, K(θ), functions. These functions are usually obtained by fitting experimental data to the corresponding θ(ψ) and K(θ) functions. The drawback of deriving the hydraulic parameters by inverse modelling is that they suffer from equifinality or non-uniqueness, and the optimal hydraulic parameters are non-physical (Pollacco et al., 2008). To reduce the non-uniqueness, it is necessary to invert the hydraulic parameters simultaneously from observations of both θ(ψ) and K(θ), and ensure the measurements cover the full range of θ from fully saturated to oven dry, which requires expensive, labour-intensive measurements.  

We present a novel procedure to derive a unique, physical set of bimodal or dual permeabilityKosugi hydraulic functions, θ(ψ) and K(θ), from inverse modelling. The Kosugi model was chosen given its parameters have direct physical meaning to the soil pore-size distribution. The challenge of using bimodal functions is they require double the number of parameters (Pollacco et al., 2017), exacerbating the problem of non-uniqueness. To address this shortcoming, we (1) derive residual soil water content from the matrix Kosugi standard deviation, (2) derive macropore hydraulic parameters from the soil water pressure boundary between macropore and matrix, and (3) dynamically constraint the matrix Kosugi hydraulic parameters. We successfully reduce the number of hydraulic parameters to optimize and constrain the hydraulic parameters without compromising the fit of the θ(ψ) and K(θ) functions.

The robustness of the methodology is demonstrated by deriving the hydraulic parameters exclusively from θ(ψ) and Ksdata, enabling satisfactory prediction of K(θ) without having measured K(θ) data. Moreover, having a reduced number of hydraulic parameters that are physical allows an improved characterization of hydraulic properties of soils prone to preferential flow, which is a fundamental issue regarding the understanding of hydrological processes.



Pollacco, J.A.P., Ugalde, J.M.S., Angulo-Jaramillo, R., Braud, I., Saugier, B., 2008. A linking test to reduce the number of hydraulic parameters necessary to simulate groundwater recharge in unsaturated soils. Adv Water Resour 31, 355–369.

Pollacco, J.A.P., Webb, T., McNeill, S., Hu, W., Carrick, S., Hewitt, A., Lilburne, L., 2017. Saturated hydraulic conductivity model computed from bimodal water retention curves for a range of New Zealand soils. Hydrol. Earth Syst. Sci. 21, 2725–2737.

How to cite: Fernández-Gálvez, J., Pollacco, J., McNeill, S., Carrick, S., Lilburne, L., Lassabatere, L., and Angulo-Jaramillo, R.: Reducing non-uniqueness of inverting bimodal soil Kosugi hydraulic parameters, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8535,, 2021.

The Daisy agroecological system model: A physically based model for preferential water flow and solute transport in drained agricultural fields
Efstathios Diamantopoulos, Maja Holbak, and Per Abrahamsen
Martin Lanzendörfer

Following the capillary bundle concept, i.e. idealizing the flow in a saturated porous media in a given direction as the Hagen-Poiseuille flow through a number of tubular capillaries, one can very easily solve what we would call the forward problem: Given the number and geometry of the capillaries (in particular, given the pore size distribution), the rheology of the fluid and the hydraulic gradient, to determine the resulting flux. With a Newtonian fluid, the flux would follow the linear Darcy law and the porous media would then be represented by one constant only (the permeability), while materials with very different pore size distributions can have identical permeability. With a non-Newtonian fluid, however, the flux resulting from the forward problem (while still easy to solve) depends in a more complicated nonlinear way upon the pore sizes. This has allowed researchers to try to solve the much more complicated inverse problem: Given the fluxes corresponding to a set of non-Newtonian rheologies and/or hydraulic gradients, to identify the geometry of the capillaries (say, the effective pore size distribution).

The potential applications are many. However, the inverse problem is, as they usually are, much more complicated. We will try to comment on some of the challenges that hinder our way forward. Some sets of experimental data may not reveal any information about the pore sizes. Some data may lead to numerically ill-posed problems. Different effective pore size distributions correspond to the same data set. Some resulting pore sizes may be misleading. We do not know how the measurement error affects the inverse problem results. How to plan an optimal set of experiments? Not speaking about the important question, how are the observed effective pore sizes related to other notions of pore size distribution.

All of the above issues can be addressed (at least initially) with artificial data, obtained e.g. by solving the forward problem numerically or by computing the flow through other idealized pore geometries. Apart from illustrating the above issues, we focus on two distinct aspects of the inverse problem, that should be regarded separately. First: given the forward problem with N distinct pore sizes, how do different algorithms and/or different sets of experiments perform in identifying them? Second: given the forward problem with a smooth continuous pore size distribution (or, with the number of pore sizes greater than N), how should an optimal representation by N effective pore sizes be defined, regardless of the method necessary to find them?

How to cite: Lanzendörfer, M.: Two sides of the same inverse problem, in identifying the pore size distribution based on experiments with non-Newtonian fluids, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14586,, 2021.

Salma Tafasca, Agnès Ducharne, and Christian Valentin

Tropical clay Oxisols have a strong granular structure, as opposed to other clay-textured soils, such as Vertisols. Their highly aggregated structure favors infiltration and drainage, while Vertisols favor runoff, due to their high content of swelling clays. Most global scale pedotransfer functions (PTFs) are derived from soils of temperate regions where Oxisols are absent, which does not allow to represent this type of tropical soils.

In ORCHIDEE, which is a state-of-the-art land surface model (LSM), soil hydraulic properties are derived from a soil texture map, using PTFs of temperate regions (Carsel and Parrish, 1988), which does not allow to represent tropical clay Oxisols. This has been shown to induce large negative evapotranspiration biases in areas effectively covered by clay Oxisols, and mapped as clay (Tafasca et al., 2020). In order to distinguish the two types of clays of different behavior, we introduce a new soil class to represent clay Oxisols. This is equivalent to splitting the clay soil texture in two sub-classes: clay Oxisols and swelling clay. First, we modify the Reynolds soil texture map (Reynolds et al., 2000) such as to represent clay Oxisols based on the FAO Soil Order Map of the World. Then we obtain the corresponding Van Genuchten parameters based on previous studies from the litterature. We evaluate the new PTFs and the modified soil texture map using the ORCHIDEE LSM.

Using the new soil texture mapping increases evapotranspiration by 11%, allowing to correct important negative bias in tropical areas. In these zones, the mean annual river discharge decreases, consistently with the observations.


How to cite: Tafasca, S., Ducharne, A., and Valentin, C.: Accounting for soil structure in pedo-transfer functions: swelling vs non swelling clays, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-10598,, 2021.

Asra Asry, Jérémie Bonneau, Gersende Fernandes, Gislain Lipeme Kouyi, Bernard Chocat, Tim D. Fletcher, and Laurent Lassabatere

Bioretention systems are increasingly used worldwide to mitigate the impacts of urban stormwater runoff on the water cycle. The proper management of bioretention systems requires accurate modeling of physical processes occurring within these systems. This study developed and tested a generic and physically-based model called Infiltron-mod. This model makes use of the Darcian approach (assuming Mualem-van Genuchten model for the description of the soil hydraulic properties) and mass conservation. The first version of the model considers evapotranspiration, overflow, exfiltration to surrounding soils, along with the filter hydraulic head and underdrain discharge. The proposed model was tested against field data from a monitored bioretention basin in Melbourne, Australia. We used two rainfall events to calibrate the model and 20 rainfall events for its validation. We achieved quite nice fits of experimental data with median NSE values in the order of 0.7-0.75 for the outflow rates. Despite good performance for outflow rates, we noticed the potential for improvement for the simulation of the height of water in the systems. Such discrepancy is probably the result of preferential flows.

As a second step, we developed a specific module to implement the dual permeability approach to model preferential flow. Such an approach may simulate the concomitancy of matrix flow in part of the system and rapid preferential infiltration into macropores. The new module Infiltron-mod-pref was implemented and investigated. Prior to its use for field data, we validated the new module against more straightforward water infiltration experiments. Several large ring infiltration tests were performed on a field dedicated to infiltrating stormwater, and the two versions of the proposed model, Infiltron-mod and Infiltron-mod-Pref. We clearly showed the benefit to account for the preferential flow in the model. The next step will be the use of Infiltron-mod_Pref for field data from the monitored bioretention basin in Melbourne.

The proposed approach then seems a useful first step to assess both performance and impact of bioretention basins for catchment-scale flow regime management and has real potential for application where user-friendly and simple model calibration and deployment are desired.

How to cite: Asry, A., Bonneau, J., Fernandes, G., Lipeme Kouyi, G., Chocat, B., Fletcher, T. D., and Lassabatere, L.: Modelling uniform and preferential flow in bioretention systems, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-10576,, 2021.

Meet the authors in their breakout text chats