The Experiments by Henry Darcy

Henry Philibert Gaspard Darcy (1803–1858) was an engineer who built the drinking water supply system of the French city of Dijon in the mid-19th century (Freeze, 1994; Brown, 2002), for which the city was deeply grateful (Philip, 1995). In doing so he developed an interest from flow of water through sands, and together with Charles Ritter experimented (in a hospital, for unclear reasons) with water flow in vertical cylinders filled with different sands “to determine the laws of flow of water through sand.” (quotation from a translated segment of Darcy (1856) in Freeze (1994)). The results were published in 1856 in an appendix to Darcy’s report on his work on Dijon’s water supply (Darcy, 1856). Darcy and Ritter installed mercury manometers at the bottom and near the top of their cylinders, and observed that the water flux density through the sand was proportional to the difference between the mercury levels. After converting mercury levels to equivalent water levels and recasting in differential form, this relationship became known as Darcy’s law, and until this day is the cornerstone of the theory of water flow in granular media like sand and silt. In such media, the voids between the grains are the pores through which the flow takes place.

Darcy identified two key properties of porous media:

This contribution traces the development of the theory of groundwater flow that emerged from Darcy’s work. It also addresses the crucial extension of the theory to water flow in soils, where part of the pore space is filled with air (“unsaturated”). The relevance of these theories for saturated and unsaturated flows is illustrated by important practical applications.

The Potential Energy of Water

The proportionality between the difference in manometer readings and the flow rate holds for laminar flow only, which occurs if the porous material is fine enough (sand, silt, clay, or a mixture of those) and the flow velocity is low. In gravels, the pores are large and flow velocities can become high enough to render the flow turbulent, and Darcy’s law loses its validity when it does. In cases where it holds, it is worthwhile to consider in some detail the equivalent water level that apparently drives the flow. This water level describes the potential energy of the water at a given location and a given time within a subsurface water body (often an aquifer). In an aquifer it can be measured simply by installing a monitoring well at the desired location: a vertical hole is drilled and a pipe inserted into the hole. A limited section (typically a meter or less) near the lower end of the pipe is perforated to create a filter section through which water can move in and out freely. The water at the filter will often have a pressure that is larger than the atmospheric pressure, and consequently, water will flow into the tube and rise until the water level in the tube is so high above the filter that the layer of water above the filter will create enough pressure to prevent more water from flowing out of the aquifer into the well. At that time, the water in the aquifer and in the well are at equilibrium, and the water level in the well accurately reflects the potential energy of the water in the aquifer around the filter (Figure 1).

Open in new tab

Figure 1. Depiction of a deep aquifer that is drained by a river. The groundwater table is higher further away from the river, indicating the water must flow toward the river. A few paths the water takes (flow lines) are indicated. The vertical features indicate monitoring wells (not to scale). The shaded lower section depicts the location of the perforated section of the well tube (the well filter), through which the water can enter and leave the well. The water level in the well therefore represents the hydraulic potential at the location of the filter. If two wells near to one another, but with filters at different depths, do not have the same water level, this indicates that the flow is not strictly horizontal.

N.B.: The vertical scale is exaggerated for clarity.

Thus, one can measure the potential energy of water at the depth of this filter section by monitoring the water level in the well with respect to a fixed level (e.g., mean sea level). By drilling multiple wells and determining their water levels with respect to the same reference level, the potential energy in a subsurface water body can be mapped. This map gives a good indication of the general flow direction within the water body. If monitoring wells with filters at different depths are installed next to each other, it is even possible to see if there is vertical flow.

The water levels thus obtained (relative to a fixed reference level) are termed hydraulic potentials or heads, because they give the potential energy of water, expressed by the water level that would give the same potential energy. Water flows in the direction that most effectively reduces its potential energy, that is, in the direction of the largest negative gradient of the hydraulic potential.

Formalization of Flows in Second-Order Partial Differential Equations

With this, Darcy’s law can be viewed as a description of a potential-driven flow: the flow rate of water is directly proportional to the negative gradient of the potential energy of the water. This positions Darcy’s law in a family of more familiar laws that are perfect analogues of each other: among those are Ohm’s law, according to which an electrical current (in Ampères) is directly proportional to the negative gradient in electrical potential (in Volts), Fick’s law of diffusion, which states that the diffusive movement of a chemical substance is directly proportional to the negative gradient in its concentration, and Fourier’s law for conduction of heat, according to which the flow of heat in a solid is directly proportional to the negative temperature gradient.

The direct proportionality between the flux density and the gradient of the hydraulic potential has major mathematical consequences that explain why Darcy’s law has become as important as it is. Darcy’s law describes motion: it explains how fast and in which direction water moves. In order to derive equations that describe subsurface water flow, this equation of motion have to be combined with the law of conservation of mass (e.g., Verruijt, 1970). In the derivation of the flow equation this involves the inspection of a control volume: a small volume that expands in each of the principal directions of the flow. For three-dimensional flows, this is a cube. If the flow tubes are parallel in one direction and the flow becomes two-dimensional, it is a strip. For one-dimensional flows (often assumed in soils when water moves predominantly vertically), it is a line. In each of the principal directions, water enters on one side of the control volume and exits at the other. The difference between the two contributes to gain or loss of water within the control volume. For sufficiently small control volumes, the difference between the flow rates can be expressed as the local gradient in the flow rate multiplied by the distance between entry and exit plane.

With the mass balance terms consisting of gradients of the flow rate, and the flow rate being proportional to the gradient in the hydraulic potential, subsurface flow equations invariably have terms in them in which the gradient of the gradient of the hydraulic potential appears. The mathematical term for these is second-order derivatives, and the resulting equations are second-order partial differential equations. Such equations arise wherever a flux density of some quantity is driven by a potential gradient, which is quite often as we saw above. The preponderance of this type of equation governing a host of natural processes (electrical current, solute diffusion, heat flow, laminar liquid flow in porous media) generated considerable interest from mathematicians, and many techniques have been developed to solve second-order partial differential equations (Carslaw & Jaeger, 1959; Bruggeman, 1999).

Theoretical Foundation of Darcy’s Law

Henry Darcy derived his law in a strictly empirical sense: he measured flow rates of water and differences in hydraulic potential (even if he did not call them that), and noticed the linear relationship between the two. It is interesting to know if Darcy’s law can be placed on a more sound theoretical footing. To do so would require a full fluid-dynamical treatment of the flow in an immensely complex labyrinth of pores, which is not realistic. For simpler geometries (flows in cylindrical pipes or between parallel plates) this is feasible though. For these flows, the friction between the pipe wall or the plates and the fluid can be accounted for. Even the internal friction within the fluid, caused by the sliding of layers of fluid over one another as they move at different velocities, can be taken into account. When external pressure, gravity, and effects of velocity gradients and changes in velocity are included in the analysis as well, the full Navier-Stokes equation emerges, which is far more complicated than Darcy’s law. The terms in this equation can be examined one by one. In doing so it emerges that the slowness of flow in porous media makes some terms very small. Groundwater rarely flows more than a few meters per day. Most often it is even much slower than that, gently traveling through the small pores between the sand grains and clay particles for a few centimeters a day, or even slower. This makes all terms in the Navier-Stokes equation that involve velocity changes and velocity gradients orders of magnitude smaller than the remaining terms. In the end, it turns out that the friction forces between the pore walls and the water, and the internal friction forces within the water, dominate the flow process. It also emerges that these forces are proportional to the overall flow velocity. They cannot be determined theoretically, because the detailed pore architecture cannot be measured, and even if it could, it would be too complex to permit the required calculations. But although the exact slope of the relationship between the gradient in the hydraulic potential and the resulting flow rate remains hidden in the pore network labyrinth, the direct proportionality between them that was observed by Darcy and generations of hydrologists and students after him follows directly from the Navier-Stokes equation, firmly embedding Darcy’s law in the theoretical framework provided by fluid dynamics.

Analytical Solutions of Groundwater Flow Problems

Jules Dupuit (1804–1866), like Henry Darcy, was a civil engineer in France, who, among other things, worked on the construction of the sewage system of Paris. (Unlike Darcy, he also was a self-taught economist of note). Dupuit and Darcy cooperated for some time (Ritzi & Boboeck, 2008) and Dupuit, like Darcy, published a monograph on subsurface flow (Dupuit, 1857). Ritzi & Bobeck (2008) argue that Darcy (1856) and Dupuit (1857) establish many of the fundamentals of subsurface flow. Dupuit introduced the differential form of Darcy’s relationship in 1863 (Nimmo & Landa, 2005), which is the form in which it is known today. It is worth noting that in his groundwater studies, Dupuit was among the first to subscribe to the idea that most of the water in an aquifer flowed slowly in a porous matrix, whereas many at the time believed that large, cavernous conduits would channel most of the flow (Brown, 2002). This opened up the way for the derivation of groundwater flow equations based on the combination of the mass conservation law with Darcy’s Law. Dupuit (1863) and Forchheimer (1886) introduced assumptions in the description of groundwater flow that allowed the resulting partial differential equation for groundwater flow to be simplified to a degree that permitted analytical solutions of it. (An analytical solution to a partial differential equation for groundwater flow is an equation that gives the value of the hydraulic potential for every point in the aquifer at every time. From this, other quantities of interest—the direction and velocity of the flow, for instance—can be easily calculated.) The Dupuit-Forchheimer assumptions are (Figure 2):

Open in new tab

Figure 2. A cross-section of the subsurface (Figure 2A), with two aquifers (composed of conductive material, such as sand and gravel) separated by a poorly conductive geological formation (termed aquitard, often composed of loam and clay). Aquifer 2 rests on a second aquitard. This could be another clay layer, or perhaps poorly fractured rock. The image depicts the flow in case the hydraulic head in aquifer 2 is larger than that of aquifer 1: this drives a small upward flow between the aquifers. For calculations, the subsurface is simplified according to the Dupuit-Forchheimer assumptions (Figure 2B). The ditches reach all the way to the first aquitard, and the variation of the aquifer thickness is so small that it is assumed to be constant. This allows the flow to be assumed horizontal.

N.B.: The vertical scale is exaggerated for clarity.

Because aquifers tend to extend over many kilometers horizontally but are mostly only tens or hundreds of meters deep, assumptions 1–3 are quite realistic in many cases. Assumption 4 is generally valid far away from a pumping well or a river that cuts into an aquifer and is fed by water flowing into it, but may not be satisfied within a few tens of meters from such a well or river.

Large aquifers are often used as sources of high-quality drinking water. Their spatial extent guarantees sufficient recharge by water percolating through the soil to permit large amounts of water to be pumped from the aquifer without excessive loss of hydraulic head. The potential yield of the wells in such an aquifer (the largest flux that the wells can sustain) can be sufficient for entire villages or small cities. Furthermore, the water in such large systems often resides underground in an anaerobic environment for hundreds or thousands of years. This means the water is generally free of microbes. The filtering effect of the porous medium through which it flows also removes particulate matter it might have carried along when it entered the aquifer as infiltrating rainfall all that time ago. The earliest studies in groundwater hydrology were concerned with drinking water wells. Thiem (1906), de Glee (1930), and Theis (1935) were the first to develop equations (solutions to partial differential equations) for flows toward drinking water wells. Their equations are still used to estimate the potential yield of such wells. The Dupuit-Forchheimer assumptions spawned numerous other analytical solutions for ever more complicated groundwater flows. In recent years mathematical tools have developed to a level of sophistication that permits complex, three-dimensional flow problems to be tackled analytically, but generally numerical approaches are used for such problems.

Determining Aquifer Properties through Pumping Tests

Later generations used the analysis of radial flows toward pumping wells pioneered by Thiem, de Glee, and Theis to develop experimental configurations of one or more pumping wells surrounded by unpumped monitoring wells to determine crucial aquifer properties, mainly effective thickness, hydraulic conductivity, or the product of the two, the transmissivity. The latter parameter gives an indication of the ease with which an aquifer can sustain a given withdrawal of water from a well. In a typical pumping test, as these experiments are called, water is pumped from the pumping well, either constantly or intermittingly. At the same time, the drop in the water levels in the monitoring wells is measured (Figure 3). Some tests keep on going until the water levels in all wells remain constant, while others measure the speed of the drawdown as a function of the distance from the pumping well and the speed of recovery when the pumping is stopped, depending on the type of flow equation that is used and the nature of the analysis that has been selected.

Open in new tab

Figure 3. A possible setup for a pumping test to determine the hydraulic properties of an aquifer. The pumping well has a filter that runs through the entire cross-section of the aquifer. The monitoring wells also have very deep filters because the extent of the drawdown of the groundwater table is not known a priori. Initially, the groundwater level is (nearly) horizontal. After pumping, a drawdown cone develops. To ensure the shape of the cone is not affected, the pumped water must be hosed off well beyond the circumference of the cone before it is allowed to infiltrate back into the soil—often several hundreds of meters away.

N.B.: The vertical scale is exaggerated for clarity

Various equations were developed for aquifers that were isolated from above and below (confined aquifers, sandwiched between impermeable layers), or on the contrary were able to exchange water with deeper aquifers (in which case they were considered leaky). Phreatic aquifers were also considered. In a complicated sequence of geological formations, these are the aquifers nearest to the soil surface (sometimes even reaching the soil surface, visible in the landscape as swamps and lakes). The groundwater table is their upper boundary. Each of these circumstances leads to a different equation that governs the flow of water toward the well. The geohydrologist can choose the extent of geohydrological surveying, the type of pumping test, the details of its execution (duration, number and spacing of monitoring wells, etc.) and the preferred computational method best suited to the conditions encountered in the field and the available budget and resources (Kruseman & de Ridder, 1990).

Numerical Solutions and Practical Groundwater Modeling

The equations governing groundwater flow can be solved numerically for very flexible conditions. Consultancy companies routinely use topographical, geological, and hydrogeological maps and hydrogeological information to build a detailed virtual representation of the subsurface: feature such as various geological formations, rivers, ditches, and drain tubes can all be included in such models. The soil can be represented (although admittedly often in a rather crude way), different land use can be incorporated, and meteorological data of many decades can drive the model. Such models are used to optimize clean-up strategies if an aquifer is contaminated, calculate how much pumping is required to drain a dug-out construction site (and at what distance to discharge the pumped water), and carry out scenario studies to see what effect future climate and land use developments have on the reliability of an aquifer as a source of drinking water, or to investigate if a proposed policy measure (e.g., to limit the application of manure to the spring and summer months) indeed has the desired effect (cleaner and clearer surface waters). Thus, in little over a century, Darcy’s law made its way from the laboratory into practical water management.

Heterogeneous Aquifers

Aquifer clean-up projects and polluted drinking water wells increased the interest of solute transport in aquifers. It quickly became apparent that not only aquifer geometry but also the spatial variation of the porosity and the hydraulic conductivity of the porous material comprising the aquifer had a major impact on the movement of dissolved contaminants in the aquifer. This heterogeneity of subsurface geological formations has spawned a large body of scientific literature that explores the aquifer-scale ramifications for water movement and solute transport of spatial variations of the porosity and the hydraulic conductivity over distances less than a meter. Many of the theoretical approaches rely on the flow lines to be long enough for mixing processes to be fully effective, which simplifies the description of the resulting distribution of solutes in the aquifer (e.g., Dagan, 1987). When mixing is less efficient and/or the complexity of the heterogeneity cannot be captured by conventional stochastic descriptions, the problem becomes much more complicated (e.g., Berkowitz, Cortis, Dentz, & Scher, 2006; Neuman, 1990).

Unfortunately, it is very difficult to obtain undisturbed samples from large depths, and hydraulic conductivity values from disturbed samples are notoriously unreliable. These studies therefore by necessity were predominantly theoretical. They offered considerable insight into the relationship between small-scale heterogeneity and large-scale aquifer properties for narrowly circumscribed statistical properties of the heterogeneities.

Experimental verification of the validity of the limitations imposed on the nature of the heterogeneities is very difficult: sampling is nearly impossible, as mentioned above, and in addition most countries prohibit the introduction of tracers in aquifers. A rare exception is the Borden aquifer in Canada: a small (~10⁴ m³), nearly rectangular and heavily instrumented aquifer in which extensive tracer experiments have been conducted. An unprecedented dense sampling scheme allowed a thorough determination of the spatial structure of the variation of the hydraulic parameters (Sudicky, 1986), which could successfully be used to predict the spreading of solutes as they traveled through the aquifer (Woodbury & Sudicky, 1991). With the recent development of noninvasive techniques such as ground penetrating radar and electrical resistivity tomography, in tandem with the computational power required to derive the spatial distribution of hydraulic properties from data generated by a very large number of sensors, an ability to see inside an aquifer is finally developing (Müller et al., 2010).

The number of monitoring wells and other sensors required to assess the movement of water and tracers or contaminants is very much larger than what is usually feasible in practical problems. In addition, the validity of the limitations imposed on the nature of the heterogeneity of the aquifer hydraulic properties is largely unknown and the data required to determine the statistical parameters of this heterogeneity will generally be lacking for practical cases of groundwater contamination. As a consequence, the tools and methods developed and used in academic aquifer research have hardly penetrated the spheres of groundwater management authorities and consultancy firms that execute actual groundwater cleanup projects.

Drainage of Agricultural Fields

A particular branch of groundwater flow theory that had a major impact on agricultural productivity is drainage theory. Drainage systems often consist of subsurface drainage tubes that discharge into a network of ditches that transport the water away from the fields. In humid climates, drainage is used to control the groundwater level. The potentially most productive agricultural lands are mostly located in regions of relatively flat topography that are often too wet to fully realize the productive potential offered by the soils and the climate.

In semi-arid areas, productivity can be boosted by irrigation. This creates an inherent risk of soil salinization, which can be controlled by leaching out salts and ensuring the groundwater table remains deep enough. In irrigated areas, the main role of drainage therefore is salinity control.

With the heavy agricultural equipment currently in use, trafficability is also important. In particular, artificially lowered groundwater levels allow tractors on the fields longer in autumn (for harvesting and plowing), but more important is early trafficability in spring, since this adds considerable flexibility in the timing of the seedbed preparation, and therefore, for the crop rotation. A drier soil also heats up more rapidly after winter, creating a longer growing season.

The flow toward drainage tubes is radial near the tubes. Furthermore, the water that needs to be discharged enters from above, so the trajectories that the groundwater water follows toward the drains (the flow lines) start from the water table. Solutions to flow equations based on the Dupuit-Forchheimer assumptions would only consider horizontal flow. This requires any water entering from above to be instantly distributed over the entire vertical cross-section so it can start flowing horizontally immediately. The radial flow toward the drain tubes cannot be represented under the Dupuit-Forchheimer assumptions either. The requirement of strictly horizontal flow means that only flow to very deep ditches (reaching all the way to the bottom of the aquifer) can be considered. All in all, this leads to an underestimation of the resistance that the flowing water encounters on its way toward the drain tube. For aquifers that are thick relative to the drain spacing (about one order of magnitude smaller), the Dupuit-Forchheimer assumptions may therefore be too crude and lead to an overestimation of the spacing of the drain tubes (e.g., Kirkham Toksöz & van der Ploeg, 1974). A full analysis of the hydraulic potential field in the vertical cross-section of the aquifer slab between the midpoint of a drain and the symmetry axis halfway between the drain tube and the next can provide more accurate estimates. This approach was pioneered by Don Kirkham in the 1950s and 1960s (see the references in Kirkham et al., 1974). The mathematical complexity of the analysis limited its application mainly to steady flow problems, but with suitable approximations some solutions for transient result were obtained (van Schilfgaarde, 1974b).

Before Kirkham, Symen Barend Hooghoudt (1901–1953) (Raats & van der Ploeg, 2005) derived an equation (Hooghoudt, 1940) that relates the desired drain spacing to the steady design infiltration rate, the saturated hydraulic conductivity of the aquifer, the aquifer thickness, and the maximum permitted elevation of the groundwater table above the level of the drain tubes halfway between the drains. The assumptions he made were similar to the Dupuit-Forchheimer assumptions, which led to an overestimation of the drain spacing. He corrected for the lack of radial flow by reducing the aquifer thickness to an effective thickness, thus increasing the resistance to the flow everywhere in the cross-section. This extra resistance resulted in an increase of the water level halfway between the drains (Figure 4). To reduce this level to the target level, the drains would have to be closer to one another. This correction was so effective that Hooghoudt’s drainage formula in combination with the expression to calculate the effective aquifer thickness is still widely used today. For layered profiles, the equation does not work though. Ernst (1954, 1956) developed a more elaborate equation that accounted for different soil layers as well as radial flow toward the drain.

Open in new tab

Figure 4. The flow toward drainage tubes in an agricultural field, and the resulting shape of the groundwater table (Figure 4A). Figure 4B depicts the simplifications used by Hooghoudt (1940). The drains have been replaced by deep ditches, the flow is considered horizontal, and the thickness d of the aquifer is smaller than the true thickness (D) to compensate for the extra resistance caused by the vertical flow component of the infiltrating water and the radial, converging flow near the drain tubes.

N.B.: The vertical scale is exaggerated for clarity.

The analyses by Kirkham, Hooghoudt, and Ernst were all concerned with steady flows, and therefore work best in moderate climates with prolonged low-intensity rainfall. If rainfall is more concentrated in heavy showers, dynamic approaches that define a drainage criterion in terms of the magnitude of the groundwater level decline and the time since the cessation of rainfall to achieve this decline are more useful. The most well-known of this is the Dumm-Glover theory (Dumm, 1954), which relied on the Dupuit-Forchheimer assumptions. It predicted an exponentially decaying water level during dry periods. The speed of the decline was characterized by the reservoir coefficient, which was proportional to the square of the drain spacing and inversely proportional to both the thickness and the hydraulic conductivity of the aquifer. Other transient approaches exist (see van Schilfgaarde, 1974b) but have not found the widespread application of the steady-state approaches and will not be further discussed.

In drainage applications too, the use of numerical computer models changed the way drainage systems are designed. A one-dimensional model of water flow in the soil profile, coupled with a plant growth model that calculates the growth of different crops as dependent on the weather and the available water, can be augmented with an approximate model for lateral flow of water toward the drains that takes into account drain spacing and drain depth. Such a model configuration can then be fed with several decades of weather data to calculate the crop yield over this period for various crop rotations, drain depths, and drain spacings. In this way, the drainage depth and spacing that gives the highest crop yield can be determined.

Extension of Darcy’s Law for Unsaturated Soils

A short distance above the groundwater table, the soil remains saturated. Higher up, air will enter the largest pores. With increasing height above the water table, increasingly smaller pores will become air-filled. In depths where roots penetrate and take up water and close to the soil surface the soil can easily dry out to such an extent that water only occurs in films around the grains and in small pockets around the contact points of grains. All regions with part of the pore space filled with air constitute the unsaturated zone. Water flow in the unsaturated zone is more complicated than groundwater flow. Research in that area was pioneered by Edgar Buckingham (1867–1940), who introduced in 1907 a version of Darcy’s law for unsaturated soils (Buckingham, 1907), possibly without being aware of Darcy’s work (Nimmo & Landa, 2005). The proportionality constant in this version would become known as the unsaturated hydraulic conductivity.

Buckingham also introduced the concept of matric potential, which expresses the degree of attraction of the water to the soil. This variable is central in unsaturated flow theory. It indicates how much force (or suction) needs to be applied to extract an infinitesimally small amount of water from the soil. At the groundwater table, the soil is saturated and the water is at atmospheric pressure. Extracting a little water does not require additional suction, and the matric potential is zero. Above the groundwater table, the water in the soil is at subatmospheric pressure and a certain force is required to extract a small amount from the soil. In accordance with the principle that water flows from high to low potential, the potential in the unsaturated region must be smaller than that at the groundwater table, and the matric potential therefore is negative in unsaturated soils. To extract the water, an even more negative potential is required. In moderately dry soils this can be achieved by applying suction to an unglazed ceramic cup buried in the soil. But soils can easily dry out to such an extent that even a vacuum cannot extract the water. Some measurement techniques measure the vapor pressure in a closed small chamber above a small soil sample. Assuming equilibrium between the liquid water in the soil and the vapor above it, the matric potential can be calculated directly from the vapor pressure.

The atmosphere typically is undersaturated with water vapor to such a degree it is at equilibrium with a very dry soil. Plants use this large difference in potential between the soil water and the atmosphere by offering a path for water that is continuous from the roots to the leaves, and is much more conductive to liquid water than the dry soil pores near the soil surface. The potential difference draws water into the roots and carries it to the leaves, where it is released into the atmosphere. With the water, the plants take up the nutrients they need to grow.

Buckingham realized that both the unsaturated hydraulic conductivity and the water content depend on the matric potential and pioneered experimental methods to observe these relationships. Thus, the elements fundamental to modern unsaturated flow theory were there: the matric potential whose gradient together with gravity drives the flow, an adapted version of Darcy’s law relating the flux density to the potential gradient, and the dependence of both the hydraulic conductivity in Darcy’s law and the soil water content on the matric potential.

The Soil Hydraulic Property Curves

A saturated porous medium can be characterized by its porosity and saturated hydraulic conductivity. But as air replaces water when the medium dries, these characteristics are merely the end points of two curves that give soil water content and the hydraulic conductivity as a function of the matric potential (Figure 5). These curves are very non-linear: neither the water content nor the hydraulic conductivity drops at a constant rate with decreasing matric potential. The curves are so closely connected to the architecture of the pore space that an experienced soil physicist can qualitatively assess the type of grains and their size distribution from looking at these curves. Sandy soils are nearly dry at matric potentials for which a heavy clay soil is still nearly saturated, for instance. Soil water characteristics (water content as a function of matric potential) can be measured relatively well. Unfortunately, the pore space architecture creates different wetting and drying patterns depending on the exact distribution of the water at the start of a wetting or drying cycle. A consequence of this is hysteresis in the soil water characteristic: it is different for every wetting-drying loop. Only the main curves, which result from drying after complete saturation and wetting after complete drying, are reproducible. This translates into a hysteretic relationship between the hydraulic conductivity and the matric potential as well. Luckily the relationship between the hydraulic conductivity and the water content exhibits much less hysteresis. To this day, measuring unsaturated hydraulic conductivities over a wide range of water contents and matric potentials is a very challenging problem that has driven considerable methodological innovations.

Open in new tab

Figure 5. An example of a hydraulic conductivity curve. The depicted data are from a sandy soil from a coastal dune in The Netherlands. Note how rapidly the hydraulic conductivity drops: from around 300 cm/day at saturation to 0.2 cm/day when the soil is at equilibrium with a groundwater table 65 cm below. The high conductivity when the soil is wet and the rapid decline are both typical for sandy soils. A clay soil would still be saturated when the groundwater would be 70 cm deeper, but the conductivity at saturation would only be a few cm per day and then drop off more gradually, making the clay soil very poorly conductive in comparison to sand for high matric potentials (close to zero) but more conductive for lower matric potentials: at -60 cm, many clay soil will conduct water better than this sand. The noisiness of the data reflects the difficulties of measuring the unsaturated hydraulic conductivity. Because its spatial variation within a few meters is even larger this is not considered a major problem.

A Partial Differential Equation for Unsaturated Flow

In 1931, Lorenzo Adolph Richards (1904–1993) presented in his PhD thesis at Cornell University and a subsequent paper on the partial differential equation that formally combined the Darcy-Buckingham law (including the matric potential gradient as the driving force of unsaturated flow) with the law of mass conservation to describe unsaturated water flow (Richards, 1931). The derivation is similar to that for the general flow equation for saturated flow. But the presence of a second mobile phase (the gas phase) required additional assumptions: Richards assumed the gas pressure to be atmospheric at all times, which requires the gas phase to be continuous (no isolated pockets of soil air) and able to flow without resistance through the soil. Without this assumption, Richards would have been forced to develop a set of coupled differential equations for the liquid and the gas phase. Furthermore, the water content and the matric potential are assumed to be at instantaneous equilibrium, which permitted the use of the soil hydraulic property curves without additional dynamical terms. Both assumptions are still widely used today. Philip (1974) elucidated the relationship between the work of Buckingham and Richards, and the brilliance of both.

Analytical Solutions of the Unsaturated Flow Equation

Richards’ equation provided a basis for systematic modeling of unsaturated flows. The extreme non-linearity of the soil water characteristic and especially the relationship between hydraulic conductivity and matric potential (hydraulic conductivity curve) proved a formidable obstacle. With computational power insufficient to deal with the non-linearities, analytical solutions were used, but these were only possible for very restricted cases. To arrive at an analytically tractable version of Richards’ equation, Klute (1952) introduced the Kirchhoff transformation in soil physics. Transformations of partial differential equations are sophisticated mathematical techniques that cast a differential equation in a particular variable into a simpler differential equation in a transformed variable. This simpler version can then be solved, and if the solution can be transformed back to the original variable (which is not always possible), one has obtained the solution to the original differential equation. In combination with the exponential hydraulic conductivity curve by Gardner (1958), the Kirchhoff transform linearized the steady-state version of Richards’ equation, which was subsequently solved for a wide variety of problems (see the review by Pullan, 1990). The exponential conductivity curve could not fully capture the full range of the hydraulic conductivity between saturation and dryness, and the limitation to steady state flow (with few exceptions requiring additional limiting assumptions) further reduced the range of applications. These limitations proved so severe that most analytical solutions have become obsolete with the advance of numerical solution methods that could more realistically model real-world conditions.

Numerical Models of Soil Hydrological Processes

Not only is soil water flow physically and mathematically much more complicated than the saturated flows in the groundwater domain below the soil, there are many more processes that interact with each other. Rainfall and snowmelt provide water; evaporation and transpiration take it away. Plants grow in soils, and their growth rate is closely connected to water availability. This creates complicated feedbacks between water supply, weather conditions, developments of the root network, which determines the capability for future water uptake, and development of the plant canopy, which strongly influences the transpiration rate.

Soil temperature becomes particularly relevant in cold climates, where frost can immobilize water and alter the structure of the pore network. Solute transport is also important, as nutrient fluxes influence growth rates and thus create another feedback through the effect of the plants on the water movement. Furthermore, the soil, unlike the groundwater, is an aerobic environment most of the time. A wide variety of chemical reactions can occur, and microbes can break down or take up many substances.

Biological activity (growing roots, burrowing animals) can change the soil structure in ways that affect the movement of water and air. The weather can do the same: frost can break large clods and crumble the soil; raindrops can destroy the clods near the surface and create a sealing crust that promotes erosion and hampers seedlings underneath from reaching the open air. Farmers purposefully change the structure of the pore network by plowing and cultivating, add organic matter through manure, and supply water, nutrients, and pesticides as needed.

Given this degree of complexity, numerical models have become very important for modeling soil processes, and are actively being used to optimize farm practices. For the models to work, the non-linearity of Richards’ equation required sophisticated and computationally expensive techniques, available since the 1980s.

These models can handle transient flows, allow a much wider range of initial and boundary conditions than analytical solutions, and can be combined with crop growth models to study the interaction between soil moisture conditions and plant growth through root water uptake. These models are widely used to optimize irrigation practices, both in scenario studies helping to determine profitable crop rotations under given climatic conditions and irrigation water availability, and in real-time use by farmers in order to optimize the timing and magnitude of individual irrigations. They also replaced to a considerable degree drainage theories based on analytical solutions of saturated flows toward drains in the design of tube or ditch drainage systems in agricultural fields in areas with high groundwater levels.

Soil Salinization and Pollution

The intensification of agriculture in Europe and North America led to imports of fodder to feed the increasing livestock of meat-producing animals. In the 1980s it became apparent that this amounted to an international transfer of nutrients. In the exporting countries this often resulted in loss of soil quality. The importing countries suffered an excess of nutrients: the manure produced was generated with fodder growing on hectares that were not part of the farms on which the animals were kept. The manure was applied on the much smaller farm area, leading to severe over-fertilizing. This was one of the causes of acid rain in western Europe (through chemical atmospheric reactions when the nitrogen escaped to the air as ammonia), and also caused widespread pollution of groundwater and surface water with nutrients.

At about the same time, the heavy use of pesticides in the decades before led to pesticides appearing in drinking water wells (Leistra & Boesten, 1989). More recently, problems arose with pharmaceuticals and hormones that are used in cattle breeding. For instance, soil bacteria are developing a resistance against antibiotics that enter the soil through manure (Udikovic-Kolic, Wichmann, Broderick, & Handelsman, 2014).

In a very different development, irrigated agriculture tends to lead to soil salinization. The irrigation water will always contain some salts. The plants only take up the water, and the salt remains in the soil. Over the years, the salt concentration in the root zone builds up, and in the long run renders the soil unproductive. In the worst cases, only one or two generations can practice agriculture, after which the soil is unproductive for many decades, if not centuries. This can be remedied to some degree by applying excess water to leach out the soils, but this is often not done adequately: if water is rare, farmers prefer to use it to produce crops. Overirrigating in order to lose some water to drainage can easily be viewed as wasting water instead of protecting the soil.

Numerous areas have saline groundwater. If these areas are irrigated, salt can enter the root zone from the top (through the irrigation water) and from below, when groundwater levels rise when excess irrigation water recharges the groundwater. In such cases, there seems no way out: irrigating as much as the crops require leads to salinization from the salts in the irrigation water. But leaching excess salts from the soil profile by periodic over-irrigation may lead to rising groundwater tables. As soon as the groundwater table reaches high enough for capillary rise to replenish water lost from the soil by evaporation from the bare soil and transpiration through the leaves of the vegetation or the crop, the root zone will be salinized from below. This process can be at least as fast as salinization from irrigation water only, and the effects can be equally disastrous. For this reason, it is imperative that an adequate drainage system (often consisting of a combination of drain tubes and ditches) is installed in such regions before crops are irrigated (e.g., van Schilfgaarde, 1974a). By discharging the water above the target level, the drainage infrastructure keeps the groundwater below the critical depth by capturing the excess irrigation water and carrying it off, away from the fields.

In the case of overfertilizing, the problem is caused by nutrients leaving the soil, whereas in the case of irrigation and pharmaceuticals/hormones in manure, the problem is caused by solutes staying in the soil. Despite their opposite nature, they put the issue of transport of solutes in soils in sharp focus. Water is the main carrier of solutes, so solute transport problems are intricately connected to the water flow. By the early 1980s it had become apparent that taking a few bags of soil into the laboratory, filling a set of columns with soil material, and then studying the movement of solutes in these columns was not providing any information that could be used for field conditions. It transpired that soils were very heterogeneous, and repacked columns would never be able to capture that.

Soil Heterogeneity

In aquifers, the porosity and the hydraulic conductivity vary in space. In soils, this spatial variation extends to the full soil water characteristic and soil hydraulic curves. Furthermore, soils with a high clay content shrink when they dry, creating a network of cracks. Soils also have large pores created by burrowing animals and decayed plant roots. At this time, this variability has not yet been fully characterized for any soil, but its effect on solute transport is well documented (e.g., Roth, Jury, Flühler, & Attinger, 1991). Those studies that tried to capture it by analyzing a large number of samples within varying soil volumes (Warrick, Mullen, & Nielsen, 1977; Hills, Hudson, & Wierenga, 1992; de Rooij, Kasteel, Papritz, & Flühler, 2004) showed an overwhelming degree of soil heterogeneity at any scale (Nielsen & Wendroth, 2003), even if cracks and biopores were not considered. Three-dimensional modeling has only taken off in the last few years, and in those, the description of soil heterogeneity was highly simplified. At this time, realistic modeling of detailed water flow in heterogeneous soils is not yet possible at the field scale.

A particular property of soils is their spatial extent: by their nature they are thin slabs that are orders of magnitude larger in the horizontal directions than in the vertical direction (Figure 6). The flow direction is predominantly vertical when the soil is unsaturated. This configuration makes it very difficult to transfer many of the stochastic approaches that have been used in groundwater hydrology to water flow and solute transport in soils. In soils, travel distances are so short that the requirements for the simplified descriptions of solute leaching that are useful for groundwater flow will hardly ever be met.

Open in new tab

Figure 6. An example of a very thin soil cover overlying ancient rock in the Hardangervidda plateau in Norway. The plateau has an area of over 6000 km², yet the soil is less than a meter thick in many places.

Various mechanisms for solute transport by flowing water have been proposed, with solute dispersion behaving similar to diffusion, entirely caused by travel time variations, or some intermediate form between these two extremes (Jury and Roth, 1990; Flühler, Durner, & Flury, 1996; Pachepsky, Benson, & Rawls, 2000). The mechanism may change as the travel distance increases. It is generally difficult to distinguish between the various mechanisms: Butters, Jury, and Ernst (1989) took soil samples until 25 m depth and de Rooij, Cirpka, Stagnitti, Vuurens, and Boll (2006) measured solute leaching at 100 locations. Nonetheless, if sufficient data of high quality are obtained in elaborate field experiments, a sophisticated data analysis can characterize the degree of solute spreading with considerable precision, even in shallow soils (Vanderborght et al., 2001; de Rooij et al., 2006).

Wetting Front Instability and Preferential Flow

Even in uniform porous media, infiltrating water can move downward in an extremely non-uniform pattern. The zone between the initially dry soil and the wetted soil (the wetting front) is said to be unstable in such cases. The wetting front will break up into distinct flow paths (often called fingers) in which the water moves down rapidly while everywhere else the wetting front will become stagnant. The flow thus bypasses a large portion of the soil, even in the absence of highly conductive cracks, wormholes, animal burrows, or channels left behind by decayed roots. Less water and nutrients carried by it are available to plants, and reactive contaminants remain in the aerobic zone shorter than they would have had the entire soil participated in the flow. Contact with organic matter (to which many contaminants will adsorb, thereby limiting their movement) is also reduced when the flow occurs through preferential flow paths. In combination, these effects increase groundwater recharge but are detrimental to crop yields and groundwater quality.

Possible causes of wetting front instability are (de Rooij, 2000)

Whatever the cause, preferential flow paths tend to develop when there is a certain resistance to overcome before infiltrating water can penetrate deeper into the soil: the matric potential needs to exceed a water-entry value before the wetting front can move further. If the hydraulic conductivity at that water-entry value is larger than the infiltration rate, a limited portion of the soil suffices to conduct the water and the wetting front can break up in distinct flow paths that do not occupy the entire soil (Raats, 1973; Hillel & Baker, 1988). In case of water repellency and air entrapment, the water-entry value can be such that the water must be at a pressure higher than the atmospheric pressure. As a result, infiltrating water tends to build up directly above the wetting front until it has built up enough pressure to penetrate further. The tips of the preferential flow paths therefore are often wetter than the section higher up. This is sometimes taken as evidence that these flows do not obey Richards’ equation. However, when the soil water retention curve has a distinct water-entry value, there is no fundamental reason why Richards’ equation should not be able to reproduce this behavior. The numerical difficulties to do so are considerable though.

If the soil is not water-repellent, low but non-zero initial water contents suffice to stabilize the wetting front (Diment & Watson, 1985) by lowering or eliminating the water-entry value. In soils that are potentially water repellent, a critical water content may exist below which a soil is actually water-repellent (Dekker, Doerr, Oostindie, Ziogas, & Ritsema, 2001). When the soil becomes wetter, it loses its water repellency and will readily take up more water. The soil is unlikely to exhibit this critical water content since it will either dry out or wet up further. If an agricultural field exhibits both low and relatively high near-surface water contents with intermediate values occurring less often, this may indicate the existence of such a critical water content, but with random variations in its value across the field. Developing water-efficient irrigation strategies is very difficult under such circumstances (Thwaites et al., 2006).

When a developing preferential flow path enters the moist zone above the groundwater table, there is no water-entry value for the water to overcome. The preferential flow path will then widen easily (Starr, Parlange, & Frink, 1986, Cho, de Rooij, & Inoue, 2005), and the negative effects of preferential flow will be less severe.

Irrigation

The salinity problems associated with irrigation were already discussed. The effect of soil salinization is so severe that none of the modern irrigation practices is unambiguously sustainable. In fact, the only historic example of sustainable irrigation is ancient Egypt, where the lands near the Nile were irrigated for thousands of years. The yearly floods by the Nile and the annual cycle of the rising and falling groundwater level as dictated by the river level created an efficient mechanism to leach salts from the soil. With the construction of dams and the introduction of measures to regulate the river level, this mechanism is no longer operative, and soil degradation now is a major problem in Egypt (Hillel, 1998, pp. 267–268). The importance of irrigation for agricultural production, and the risks for future generation associated with soil salinization are well recognized: several countries have research institutes devoted to it. Some of those are the U.S. Salinity Laboratory in California (of which van Schilfgaarde was a director, and where Richards worked), the Volcani Center in Israel, and the International Water Management Institute in Sri Lanka.

There are two distinct developments in irrigation: the development of highly efficient methods of applying water to the soil, and the reuse of water for irrigation. Early irrigation by furrows that were periodically flooded distributed water unequally over the field. Sprinklers emitted droplets that traveled through the air over many meters and then were partially intercepted by the vegetation, which caused considerable water loss through evaporation and wind drift. Modern techniques have a more precise application of water. Pivot irrigation systems still use sprinklers, but they move in circles around a central pivot, with the nozzles close to the crop. This allows for a more precise dosage and reduces losses. Even more efficient are microirrigation systems in which water and fertilizer are delivered drop by drop directly in the root zone. Daniel Hillel (b. 1930), actively researching in both Israel and the United States, was one of the main proponents of this technique and was awarded the World Food Prize for this in 2012.

The reuse of water for irrigation nowadays focuses mainly on using raw or treated municipal sewage water for irrigation. For crops like lettuce, the leaves of which are digested directly, microbial contamination can be a risk factor. For other crops, one problem that has emerged is that prolonged irrigation with treated sewage water can render the soil hydrophobic (Wallach, Ben-Arie, & Graber, 2005). As a consequence, water will not be absorbed as easily, which obviously limits the effectiveness of irrigation. Furthermore, water repellency can lead to the formation of preferential flow paths.

Groundwater Recharge

Groundwater recharge is net flow from the soil to the groundwater. It sustains the groundwater resources that are often tapped to provide drinking water. If a river cuts through an aquifer, groundwater recharge is in part discharged into these streams, thereby securing a certain flow rate, which has important effects on river ecology, among other things. Even in some relatively humid areas, groundwater recharge is diminishing. Modern agriculture is suspected to play a role. Crop yields are more or less proportional to the total amount of water transpired by the crop in its growth cycle. More intensive use of fertilizers and more efficient pest control therefore tend to increase transpiration, leaving less water available for recharging the groundwater. In forestry, planting evergreen coniferous trees, which transpire year-round, also increases transpiration. Simultaneously, ever increasing amounts of groundwater are pumped up for domestic and industrial use. In combination, this can lead to groundwater droughts, even in countries with plenty of rainfall.

The problem is more pressing in hot semi-arid regions. With increasing population density, water scarcity is an ever increasing problem there. The communities in these regions are often supported by freshwater from aquifers, because rainfall is irregular and rivers are ephemeral. If such communities are (semi-)nomadic or practice sustenance farming they tend to have no disposable income that allows them to buy whatever water and food they need, and can only sustain themselves if they do not deplete the aquifer that supports them.

Unfortunately it has proven extremely difficult to estimate groundwater recharge of such aquifers (Scanlon et al., 2006). They are likely to be fed by rare heavy rainfall events. Rainfall data to quantify the occurrence of such events are often lacking. It has also proven very difficult to calculate the percolation of infiltrating water toward a groundwater table that is many meters deep. The way down can take years or decades, and a considerable amount of water is likely to evaporate during that time (de Vries & Simmers, 2002). Several types of problems present themselves: Richards’ equation ignores gas flow, but in this problem, diffusion of water vapor is of crucial importance. A more complete description of liquid water flow, gas flow, vapor diffusion, and heat flow is necessary to attack the problem. In addition, soil physics has mostly been concerned with relatively moist soils, since the bulk of the water movement takes place when the soil is wet. For groundwater recharge this still applies, but because the dry periods can last for months, years, or decades, they gain importance. This means that more attention needs to be paid to the distribution of liquid water in soils at low water contents, not only for the liquid flow, but also for the vapor flow. Finally, setting up long-term reliable weather monitoring stations and water content sensor networks is extremely difficult in remote areas with a very hostile climate. It remains to be seen if models and other tools that rely on Darcy’s law can perform in this arena.

Physical Limitations

Darcy’s law (including Buckingham’s version) requires the hydraulic potential to be continuous and the porous matrix to be rigid. In saturated flow, these conditions are mostly met. But in swelling and shrinking clay soils, not only does the pore space occupied by water and air change as the soils wets up or dries out, but the pore space itself changes in volume and shape. The added complications notwithstanding, flows in such media can be modeled in a suitably expanded Darcian framework (e.g., Greco, 2002), with the cautionary note that the bulk of the flow occurs through the cracks and will generally be turbulent and therefore non-Darcian.

Discontinuity of the hydraulic potential occurs in soils that are so dry that the liquid phase is limited to thin films around the grains and in pendular rings around the contact points where grains touch (Bear & Bachmat, 1991, pp. 336–340). When the films are thin enough, the attracting forces between the solids and the water are such that the water is adsorbed to the soil matrix in varying degrees and no longer behaves like a free liquid. The matric potential and the (liquid) water content are poorly defined under such circumstances, and the Darcy-Buckingham law does not perform well. When water infiltrates in such a dry soil, the wetting process proceeds incrementally, as pores are filled in sudden Haines jumps (Hillel, 1998, pp. 160–161), in which the waster nearly instantaneously fills a void between narrow pore throats (Lu, Biggar, & Nielsen, 1994). This process is decidedly non-Darcian. Immediately behind the wetting front, the jumps are rapidly smoothed out and Richards’ equation generally describes the flow well, provided the soil hydraulic properties account for the water-entry value.

In extreme cases, the water-entry value can be so large that water must have a pressure larger than the atmospheric pressure. In those cases, the soil is water-repellent, often caused by coating of organic substances on the grains. As discussed above this hampers infiltration and can lead to preferential flow. Although non-Darcian on first appearance, this type of flow can be handled by Richards’ equation if the soil hydraulic properties include a distinct water-entry value. For soils dryer than that, the hydraulic conductivity should be set to zero in front of an advancing wetting front. Infiltrating water is forced to stagnate until the pressure build-up suffices to overcome the water-entry value. However, it cannot be expected that a numerical model with this setup will predict the size of the preferential flow paths correctly.

Scale Limitations: The Regional, Continental, and Global Water Cycle

Early applications of Darcy’s law were concerned with flow toward drinking water wells. Later, the flow toward drains and in agricultural fields became of interest. Modern-day water management requires quantification of flows on the scale of catchments. The terrestrial part of the hydrological cycle is of tremendous importance for climate change research: the global distribution of water is changing, and with it the rainfall patterns. Evaporation and transpiration have a considerable effect on atmospheric temperatures: solar energy used to convert liquid water to vapor is no longer available to be converted to heat. All of this requires studies of groundwater flows over large distances—much larger than the local flow toward a well. The hydrology of soils is currently well understood at the scale of a square meter and reasonably well at the scale of a field. But there is an urgent need to quantify it on scales of basins and continents to make meaningful contributions to climate change research and to measures aimed at maintaining biodiversity and agricultural productivity even if climate zones shift.

Darcy’s law is fundamentally based on a continuum view of a porous medium. Within the range of validity of Darcy’s law, there is no need to have a distinct view of the pore architecture, because it concerns itself with much larger scales. On the other hand, the hydraulic potential that drives the flow must be continuous, otherwise its gradient that drives the flow does not exist. Also, the flow must be proportional to that gradient. This proportionality cannot be assumed to exist over an arbitrary range of scales (de Rooij, 2009). If it breaks down at some scale, the second-order partial differential equations that have proven so valuable in describing subsurface flow in the groundwater and in the soil no longer have meaning. Evidence of this can be found in large-scale groundwater models that work with spatial discretizations of tens of meters. In such models it is observed that the hydraulic conductivity defined for a volume of 10³ to 10⁵ m³ becomes dependent upon the boundary conditions (the conditions imposed on the outer edges of the aquifer) and external forcings (such as the infiltration rate from the soil above), which is in violation of Darcian groundwater flow theory. This may be attributable to deficiencies in the discretization, but could very well be a manifestation of the breakdown of the validity of the Darcian theory at that scale. At that point we have reached the limit of Darcy’s law. We would be forced to describe the subsurface of a country or a continent with maps of soil hydraulic properties and aquifer properties with a resolution of a few meters at most for the soil, and perhaps tens of meters for aquifers, even for an approximate simulation. Endeavors to develop theories of landscape evolution that can generate such maps are being undertaken, while simultaneously a quest has begun for an alternative hydrological theory specifically aimed at these super-Darcian scales.