Unsaturated Zone Flow Processes and Aquifer Response Time in the Chalk Aquifer,Brighton, South East England

The Chalk aquifer is one of the main sources of water in South East England. The unsaturated zone in the aquifer plays an important role controlling the time and magnitude of recharge and is major pathway for contaminant transport to the water table. A range of previous work has addressed flow processes in the Chalk unsaturated zone, but physical understanding is still incomplete. Here we present the results of a study on flow mechanism in the Chalk unsaturated zone using a combination of statistical analysis and novel laboratory methods. The study was undertaken at three sites (North Heath Barn [NHB], Pyecombe East [PE], and Preston Park [PP]) on the Chalk of the Brighton block, South East England. Daily and hourly time series data of groundwater level and rainfall were correlated. The results show that a slower groundwater level response to rainfall occurs during dry seasons (summer and autumn) when the amount of effective rainfall is less than 4 mm/day, with a thicker and drier unsaturated zone. A faster response occurs during wet seasons (winter and spring) when the daily effective rainfall exceeds 4 mm/day with a thinner and wetter unsaturated zone. Periods of very rapid response (within 15 h) were observed during wet seasons at NHB and PE sites, with unsaturated hydraulic conductivity (K_u) inferred to reach 839 mm/day. A slower response was observed at an urbanized site (PP) as a result of reduction in direct recharge due to reduced infiltration, due to presences of impermeable infrastructure covering the area around PP borehole. Laboratory measurements of K_u of the Chalk matrix using a geotechnical centrifuge show variation from 4.27 to 0.07 mm/day, according to the level of saturation. Thus, the rapid responses cannot be linked to matrix flow only but indicate the contribution of fracture and karstic flow processes in conducting water.     

Nitrate fluctuations at the water table: implications for recharge processes and solute transport in the Chalk aquifer

This study investigates fluctuations in nitrate concentration at the water table to improve understanding of unsaturated zone processes in the Chalk aquifer. Sampling was conducted using a novel multi‐level sampler during periods of water table rise over 5 years at a vertical resolution of 0.05 m. Nitrate concentration increased as the water table seasonally recovered, with similar inter‐annual trends with depth. The rising water table activated horizontal fractures facilitating the delivery of water elevated by up to 10 mg/l of nitrate with respect to the adjacent groundwater below. These fractures are considered to activate via piston displacement of water from the adjoining matrix. Hydrograph analysis identified 16 events which significantly perturbed the water table within 24–48 h of rainfall. Consistent nitrate concentrations indicate recharge through persistent fracture flow from the surface was not generally the primary driver of the rapid water table response during these events. Instead, the response was attributed to the piston displacement of porewater immediately above the water table. However, a single event in November 2012 delivered relatively dilute recharge indicating rapid persistent fracture flow following rainfall was possible to a depth of 14–15 m. Decreases in porewater nitrate concentration around fracture horizons and the dilution of many groundwater samples with respect to porewaters indicate a fresher source of water at depth. This was considered most likely to be a result of near surface water bypassing the matrix because of widespread mineralization on fracture surfaces, which retard water and solute exchange. Therefore, persistent fracture flow maybe considered a frequent process, operating independently of the matrix, and is not necessarily event driven. Copyright © 2015 John Wiley & Sons, Ltd.     

Prediction of the arrival of peak nitrate concentrations at the water table at the regional scale in Great Britain

A simple process‐based approach to predict regional‐scale loading of nitrate at the water table was implemented in a GIS for Great Britain. This links a nitrate input function, unsaturated zone thickness, and lithologically dependent rate of nitrate unsaturated zone travel to estimate arrival time of nitrate at the water table. The nitrate input function is the loading at the base of the soil and has been validated using unsaturated zone porewater profiles. The unsaturated zone thickness uses groundwater levels based on regional‐scale observations infilled by interpolated river base levels. Estimates of the rate of unsaturated zone travel are attributed from regional‐scale hydrogeological mapping. The results indicate that peak nitrate loading may have already arrived at the water table for many aquifers, but that it has not where the unsaturated zone is relatively thick There are contrasting outcomes for the two main aquifers which have similar unsaturated zone velocities, the predominantly low relief Permo‐Triassic sandstones, and the Chalk, which forms significant topographic features. For about 60% of the Chalk, the peak input has not yet reached the water table and will continue to arrive over the next 60 years. The methodology is readily transferable and provides a robust method for estimating peak arrival time for any diffuse conservative pollutant where an input function can be defined at a regional scale and requires only depth to groundwater and a hydrogeological classification. The methodology is extendable in that if additional information is available this can easily be incorporated into the model scheme. British Geology Survey © NERC 2011. Hydrological Process © 2011 John Wiley & Sons, Ltd     

Variable-density groundwater flow and solute transport in irregular 2D fracture networks

Numerical simulations of variable-density flow and solute transport have been conducted to investigate dense plume migration for various configurations of 2D fracture networks. For orthogonal fractures, simulations demonstrate that dispersive mixing in fractures with small aperture does not stabilize vertical plume migration in fractures with large aperture. Simulations in non-orthogonal 2D fracture networks indicate that convection cells form and that they overlap both the porous matrix and fractures. Thus, transport rates in convection cells depend on matrix and fracture flow properties. A series of simulations in statistically equivalent networks of fractures with irregular orientation show that the migration of a dense plume is highly sensitive to the geometry of the network. If fractures in a random network are connected equidistantly to the solute source, few equidistantly distributed fractures favor density-driven transport. On the other hand, numerous fractures have a stabilizing effect, especially if diffusive transport rates are high. A sensitivity analysis for a network with few equidistantly distributed fractures shows that low fracture aperture, low matrix permeability and high matrix porosity impede density-driven transport because these parameters reduce groundwater flow velocities in both the matrix and the fractures. Enhanced molecular diffusion slows down density-driven transport because it favors solute diffusion from the fractures into the low-permeability porous matrix where groundwater velocities are smaller. For the configurations tested, variable-density flow and solute transport are most sensitive to the permeability and porosity of the matrix, which are properties that can be determined more accurately than the geometry and hydraulic properties of the fracture network, which have a smaller impact on density-driven transport.     

Analytical Solutions for Water Flow and Solute Transport in the Unsaturated Zone

Analytical solutions for the water flow and solute transport equations in the unsaturated zone are presented. We use the Broadbridge and White nonlinear model to solve the Richards’ equation for vertical flow under a constant infiltration rate. Then we extend the water flow solution and develop an exact parametric solution for the advection-dispersion equation. The method of characteristics is adopted to determine the location of a solute front in the unsaturated zone. The dispersion component is incorporated into the final solution using a singular perturbation method. The formulation of the analytical solutions is simple, and a complete solution is generated without resorting to computationally demanding numerical schemes. Indeed, the simple analytical solutions can be used as tools to verify the accuracy of numerical models of water flow and solute transport. Comparison with a finite-element numerical solution indicates that a good match for the predicted water content is achieved when the mesh grid is one-fourth the capillary length scale of the porous medium. However, when numerically solving the solute transport equation at this level of discretization, numerical dispersion and spatial oscillations were significant.     

Analysis of water-level response to rainfall and implications for recharge pathways in the Chalk aquifer, SE England

Water table response to rainfall was investigated at six sites in the Upper, Middle and Lower Chalk of southern England. Daily time series of rainfall and borehole water level were cross-correlated to investigate seasonal variations in groundwater-level response times, based on periods of 3-month duration. The time lags (in days) yielding significant correlations were compared with the average unsaturated zone thickness during each 3-month period. In general, for cases when the unsaturated zone was greater than 18 m thick, the time lag for a significant water-level response increased rapidly once the depth to the water table exceeded a critical value, which varied from site to site. For shallower water tables, a linear relationship between the depth to the water table and the water-level response time was evident. The observed variations in response time can only be partially accounted for using a diffusive model for propagation through the unsaturated matrix, suggesting that some fissure flow was occurring. The majority of rapid responses were observed during the winter/spring recharge period, when the unsaturated zone is thinnest and the unsaturated zone moisture content is highest, and were more likely to occur when the rainfall intensity exceeded 5 mm/day. At some sites, a very rapid response within 24 h of rainfall was observed in addition to the longer term responses even when the unsaturated zone was up to 64 m thick. This response was generally associated with the autumn period. The results of the cross-correlation analysis provide statistical support for the presence of fissure flow and for the contribution of multiple pathways through the unsaturated zone to groundwater recharge.     

Temporal variability of solute transport under vadose zone conditions

The heterogeneity of the solute flux field in the horizontal plane at the field scale has been documented in several field studies. On the other hand, little information is available on the persistence of certain solute transport scenarios over consecutive infiltration cycles. This study was initiated to analyse the recurrence of solute leaching behaviour as estimated in two soil column tests emphasizing the preferential flow phenomenon. Twenty-four small-sized soil samples were subjected to two consecutive unsaturated steady-state flow leaching experiments with bromide as tracer. Observed breakthrough curves (BTCs) were analysed by the method of moments and by the advection–dispersion equation (ADE) to classify solute behaviour. Frequency distributions of the parameters indicating the solute velocity were heavily skewed or bimodal, reflecting the broad variability of the leaching scenarios, including some with pronounced preferential solute breakthrough. Exclusion of the preferential flow columns from our calculations revealed an average amount of 37% of immobile water. The large-scale BTCs derived from assembling the individual concentration courses of each run showed similar features, such as an early bromide breakthrough. However, two distinct apices, viz. one preferential and one matrix, were observed only in the first run, whereas the concentration decrease between the peaks was missing from the second run. A change in soil structure with continuous leaching was presumed to modify the interplay of the various flow domains, thereby altering the spreading of the BTCs. Correlation analysis between parameters of both tests suggests that preferential transport conditions are likely to occur at the same locations in the field over several infiltration cycles, whereas the ‘classical’ or expected matrix flow is time variant and therefore seems to be hardly predictable. © 1998 John Wiley & Sons, Ltd.     

A new lumped-parameter model for flow in unsaturated dual-porosity media

A new lumped-parameter approach to simulating unsaturated flow processes in dual-porosity media such as fractured rocks or aggregated soils is presented. Fluid flow between the fracture network and the matrix blocks is described by a non-linear equation that relates the imbibition rate of the local difference in liquid-phase pressure between the fractures and the matrix blocks. Unlike a Warren-Root-type equation, this equation is accurate in both the early and late time regimes. The fracture/matrix interflow equation has been incorporated into an existing unsaturated flow simulator, to serve as a source/sink term for fracture gridblocks. Flow processes are then simulated using only fracture gridblocks in the computational grid. This new lumped-parameter approach has been tested on two problems involving transient flow in fractured/porous media, and compared with simulations performed using explicit discretisation of the matrix blocks. The new procedure seems to accurately simulate flow processes in unsaturated fractured rocks, and typically requires an order of magnitude less computational time than do simulations using fully-discretised matrix blocks.     

A quasi three-dimensional model for flow and transport in unsaturated and saturated zones: 1. Implementation of the quasi two-dimensional case

A quasi three-dimensional (QUASI 3-D) model is presented for simulating the subsurface water flow and solute transport in the unsaturated and in the saturated zones of soil. The model is based on the assumptions of vertical flow in the unsaturated zone and essentially horizontal groundwater flow. The 1-D Richards equation for the unsaturated zone is coupled at the phreatic surface with the 2-D flow equation for the saturated zone. The latter was obtained by averaging 3-D flow equation in the saturated zone over the aquifer thickness. Unlike the Boussinesq equation for a leaky-phreatic aquifer, the developed model does not contain a storage term with specific yield and a source term for natural replenishment. Instead it includes a water flux term at the phreatic surface through which the Richards equation is linked with the groundwater flow equation. The vertical water flux in the saturated zone is evaluated on the basis of the fluid mass balance equation while the horizontal fluxes, in that equation, are prescribed by Darcy law. A 3-D transport equation is used to simulate the solute migration. A numerical algorithm to solve the problem for the general quasi 3-D case was developed. The developed methodology was exemplified for the quasi 2-D cross-sectional case (QUASI2D). Simulations for three synthetic problems demonstrate good agreement between the results obtained by QUASI2D and two fully 2-D flow and transport codes (SUTRA and 2DSOIL). Yet, simulations with the QUASI2D code were several times faster than those by the SUTRA and the 2DSOIL codes.     

Eulerian spatial moments for solute transport in three-dimensional heterogeneous,dual-permeability media

A Eulerian analytical method is developed for nonreactive solute transport in heterogeneous, dual-permeability media where the hydraulic conductivities in fracture and matrix domains are both assumed to be stochastic processes. The analytical solution for the mean concentration is given explicitly in Fourier and Laplace transforms. Instead of using the fast fourier transform method to numerically invert the solution to real space (Hu et al., 2002), we apply the general relationship between spatial moments and concentration (Naff, 1990; Hu et al., 1997) to obtain the analytical solutions for the spatial moments up to the second for a pulse input of the solute. Owing to its accuracy and efficiency, the analytical method can be used to check the semi-analytical and Monte Carlo numerical methods before they are applied to more complicated studies. The analytical method can be also used during screening studies to identify the most significant transport parameters for further analysis. In this study, the analytical results have been compared with those obtained from the semi-analytical method (Hu et al., 2002) and the comparison shows that the semi-analytical method is robust. It is clearly shown from the analytical solution that the three factors, local dispersion, conductivity variation in each domain and velocity convection flow difference in the two domains, play different roles on the solute plume spreading in longitudinal and transverse directions. The calculation results also indicate that when the log-conductivity variance in matrix is 10 times less than its counterpart in fractures, it will hardly influence the solute transport, whether the conductivity field is matrix is treated as a homogeneous or random field.     

A test case for the simulation of three-dimensional variable-density flow and solute transport in discretely-fractured porous media

A test case has been developed for three-dimensional simulations of variable-density flow and solute transport in discretely-fractured porous media. The simulation domain is a low-permeability porous matrix cube containing a single non-planar fracture. The initial solute concentration is zero everywhere. A constant solute concentration is assigned to the top of the domain, which increases near-top fluid density and induces downward density-driven flow. The test case is therefore comparable to downwelling of a dense brine below a saline disposal basin or a waste repository. Numerous fingers and distinct convection cells develop early in the fracture but the fingers later coalesce and convection becomes less apparent. To help test other variable-density flow and transport models, results of the test case are presented both qualitatively (concentration contours and velocity fields) and quantitatively (penetration depth, mass flux, total mass stored, maximum fracture and matrix velocity).     

Evaluating Hydraulic Parameters of Flow and Transport of a Nonreactive Solute in a River Regulated by a Chain of Hydropower Plants

A simple and efficient technique is proposed for use in combination with a known method for solving the Saint-Venant equation and the equation of advection–diffusion transport of a solute. This technique allows an overall calculation when there are several hydropower plant dams in a river channel. Practical calculations have been made for the upper reach of the Vuoksa River for the most realistic scenarios. The results of the calculations supplement the scarce field data on water level elevations, water flow, and mean flow velocities for quasi-steady-state conditions and for cases of water releases spreading downriver. The effect of water releases on the distribution of solute concentration along the river is assessed and the dependence of the pattern of a conventional pollutant spreading along the river channel after an impulse discharge on the coefficient of longitudinal dispersion is examined.     

Changes in water table level influence solute transport in uniform porous media

Changes in the water table level result in variable water saturation and variable hydrological fluxes at the interface between the unsaturated and saturated zone. This may influence the transport and fate of contaminants in the subsurface. The objective of this study was to examine the impact of a decreasing and an increasing water table on solute transport. We conducted tracer experiments at downward flow conditions in laboratory columns filled with two different uniform porous media under static and transient flow conditions either increasing or decreasing the water table. Tracer breakthrough curves were simulated using a mobile–immobile transport model. The resulting transport parameters were compared to identify dominant transport processes. Changes in the water table level affected dispersivities and mobile water fractions depending on the direction of water table movement and the grain size of the porous media. In fine glass beads, the water flow velocity was similar to the decline rate of the water table, and the mobile water fraction was decreased compared with steady‐state saturated conditions. However, immobile water was negligible. In coarse glass beads, water flow was faster because of fingered flow in the unsaturated part, and the mobile water fraction was smaller than in the fine material. Here, a rising water table led to an even smaller mobile water fraction and increased solute spreading because of diffusive interaction with immobile water. We conclude that changes of the water table need to be considered to correctly simulate transport in the subsurface at the transition of the unsaturated–saturated zone. Copyright © 2014 John Wiley & Sons, Ltd.     

Catchment‐scale conceptual modelling of water and solute transport in the dual flow system of the karst critical zone

Hydrological and hydrochemical processes in the critical zone of karst environments are controlled by the fracture‐conduit network. Modelling hydrological and hydrochemical dynamics in such heterogeneous hydrogeological settings remains a research challenge. In this study, water and solute transport in the dual flow system of the karst critical zone were investigated in a 73.5‐km² catchment in southwest China. We developed a dual reservoir conceptual run‐off model combined with an autoregressive and moving average model with algorithms to assess dissolution rates in the “fast flow” and “slow flow” systems. This model was applied to 3 catchments with typical karst critical zone architectures, to show how flow exchange between fracture and conduit networks changes in relation to catchment storage dynamics. The flux of bidirectional water and solute exchange between the fissure and conduit system increases from the headwaters to the outfall due to the large area of the developed conduits and low hydraulic gradient in the lower catchment. Rainfall amounts have a significant influence on partitioning the relative proportions of flow and solutes derived from different sources reaching the underground outlet. The effect of rainfall on catchment function is modulated by the structure of the karst critical zone (e.g., epikarst and sinkholes). Thin epikarst and well‐developed sinkholes in the headwaters divert more surface water (younger water) into the underground channel network, leading to a higher fraction of rainfall recharge into the fast flow system and total outflow. Also, the contribution of carbonate weathering to mass export is also higher in the headwaters due to the infiltration of younger water with low solute concentrations through sinkholes.     

Improvement of the solute transfer in a conceptual unsaturated zone scheme: a case study of the Seine River basin

For predicting the evolution of solute concentrations in groundwater and testing the impact of remediation policies, a coupling between the agronomical model STICS and the hydrogeological model MODCOU was implemented. When applied to the Seine River basin, this model accurately represents the temporal evolution of average nitrate concentrations in the aquifer, but with large local errors. We propose an improvement to the simple unsaturated zone (UZ) scheme NonsatSW used in STICS–MODCOU. The modifications are based on a comparison with the mechanistic model Metis considered as a reference as it solves Richards' equation. A more realistic saturation profile and a varying percolation rate are integrated in NonsatSW. This new model, named NonsatVG, is assessed by comparing it with NonsatSW and Metis. In an ideal case, NonsatVG generates a solute transfer and a dispersion closer to that of Metis than of NonsatSW. In real cases, without additional calibration, NonsatVG and Metis simulate better the average transfer velocities of the observed nitrate profiles. Furthermore, modifications in NonsatVG give a direct relationship between the depth of the water table and the saturation profile. We obtain, therefore, as in Metis, an evolution of the solute transfer velocity depending on the piezometric level. These dynamics are not simulated in NonsatSW. Despite a modified water transfer through the UZ, NonsatVG is also as valid as NonsatSW in the modelling of water transfer to the saturated zone. Finally, an application to the Seine basin shows that solute transfer velocities are lower with NonsatVG than with NonsatSW, but are in better agreement with literature. Copyright © 2010 John Wiley & Sons, Ltd.     

Modeling the Impact of Stream Discharge Events on Riparian Solute Dynamics

The biogeochemical composition of stream water and the surrounding riparian water is mainly defined by the exchange of water and solutes between the stream and the riparian zone. Short-term fluctuations in near stream hydraulic head gradients (e.g., during stream flow events) can significantly influence the extent and rate of exchange processes. In this study, we simulate exchanges between streams and their riparian zone driven by stream stage fluctuations during single stream discharge events of varying peak height and duration. Simulated results show that strong stream flow events can trigger solute mobilization in riparian soils and subsequent export to the stream. The timing and amount of solute export is linked to the shape of the discharge event. Higher peaks and increased durations significantly enhance solute export, however, peak height is found to be the dominant control for overall mass export. Mobilized solutes are transported to the stream in two stages (1) by return flow of stream water that was stored in the riparian zone during the event and (2) by vertical movement to the groundwater under gravity drainage from the unsaturated parts of the riparian zone, which lasts for significantly longer time (> 400 days) resulting in long tailing of bank outflows and solute mass outfluxes. We conclude that strong stream discharge events can mobilize and transport solutes from near stream riparian soils into the stream. The impact of short-term stream discharge variations on solute exchange may last for long times after the flow event.     

Mobility of major ions and nutrients in the unsaturated zone during paddy cultivation: a field study and solute transport modelling approach

Study of the movement of water and solute within soil profiles is important for a number of reasons. Accumulation of prominent contaminants from agricultural chemicals in the unsaturated zone over the years is a major concern in many parts of the world. As a result, the unsaturated zone has been a subject of great research interest during the past decade. Hence, an intensive field study was conducted in a part of Palar and Cheyyar river basins to understand the variation of major ions and nutrients in the soil zone during paddy cultivation. The chloride and nitrate data were used to model the movement of these chemicals in the unsaturated zone using the HYDRUS‐2D model. The field study shows that fertilizer application and irrigation return flow increases the major ions and nutrients concentration in the unsaturated zone. Further, the nutrient concentrations are regulated by plant uptake, fertilizer application and infiltration rate. Additionally, denitrification and soil mineralization processes also regulate the nitrogen concentration in the unsaturated zone. The solute transport modelling study concluded that the simulated results match reasonably with the observed trends. Simulated concentrations of chloride and nitrate for a 5‐year period indicate that the concentrations of these ions fluctuate in a cyclic manner (from 60 to 68 mg l^?1 and from 3·4 to 3·5 mg l^?1 respectively in groundwater) with no upward and downward trend. The influence of excessive fertilizer application on groundwater was also modelled. The model predicts an increase of about 17 mg l^?1 of chloride and 2·3 mg l^?1 of nitrogen in the groundwater of this area when the application of fertilizers is doubled. The model indicates that the present level of use of agrochemicals is no threat to the groundwater quality. Copyright © 2007 John Wiley & Sons, Ltd.     

Ecohydrologic considerations for modeling of stable water isotopes in a small intermittent watershed

Naturally occurring stable water isotope tracers provide useful information for hydrologic model development and calibration. Existing models include varied approaches concerning unsaturated zone percolation mixing (preferential versus matrix flow) and evapotranspiration (ET) partitioning. We assess the impact of unsaturated zone simplifying assumptions when simulating the Shale Hills Watershed, a small (7.9 ha), temperate, forested watershed near Petersburg, Pennsylvania, USA, with a relatively simple model. We found that different model structures/assumptions and parameterizations of unsaturated zone percolation had substantial impacts on the agreement between simulated and observed unsaturated‐zone water isotopic signatures. We show that unsaturated zone percolation mixing primarily affects the unsaturated zone δ¹⁸O and δ²H during winter and spring and that percolation was best represented as a combination of both preferential and matrix flow. We evaluate the importance and implications related to the partitioning of ET into evaporation and transpiration and demonstrated that incorporation of a plant growth model for ET partitioning substantially improved reproduction of observed hydrologic isotopic patterns of the unsaturated zone during the spring season. We show that unsaturated zone percolation mixing and ET partitioning approaches do not substantially influence stream δ¹⁸O and δ²H and conclude that observed streamflow isotopic data is not always a strong predictor of model performance with respect to intrawatershed processes.     

Rapid imbibition of water in fractures within unsaturated sedimentary rock

The spontaneous imbibition of water and other liquids into gas-filled fractures in variably-saturated porous media is important in a variety of engineering and geological contexts. However, surprisingly few studies have investigated this phenomenon. We present a theoretical framework for predicting the 1-dimensional movement of water into air-filled fractures within a porous medium based on early-time capillary dynamics and spreading over the rough surfaces of fracture faces. The theory permits estimation of sorptivity values for the matrix and fracture zone, as well as a dispersion parameter which quantifies the extent of spreading of the wetting front. Quantitative data on spontaneous imbibition of water in unsaturated Berea sandstone cores were acquired to evaluate the proposed model. The cores with different permeability classes ranging from 50 to 500 mD and were fractured using the Brazilian method. Spontaneous imbibition in the fractured cores was measured by dynamic neutron radiography at the Neutron Imaging Prototype Facility (beam line CG-1D, HFIR), Oak Ridge National Laboratory. Water uptake into both the matrix and the fracture zone exhibited square-root-of-time behavior. The matrix sorptivities ranged from 2.9 to 4.6 mm s^−0.5, and increased linearly as the permeability class increased. The sorptivities of the fracture zones ranged from 17.9 to 27.1 mm s^−0.5, and increased linearly with increasing fracture aperture width. The dispersion coefficients ranged from 23.7 to 66.7 mm² s⁻¹ and increased linearly with increasing fracture aperture width and damage zone width. Both theory and observations indicate that fractures can significantly increase spontaneous imbibition in unsaturated sedimentary rock by capillary action and surface spreading on rough fracture faces. Fractures also increase the dispersion of the wetting front. Further research is needed to investigate this phenomenon in other natural and engineered porous media.