Modeling of Flow and Transport in Saturated and Unsaturated Porous Media

Abstract

Modeling fluid flow and transport processes in porous media is a relevant topic for a wide range of applications. In water resources problems, this topic presents specific challenges related to the multiphysical processes, large time and space scales, heterogeneity and anisotropy of natural porous media, and complex mathematical models characterized by coupled nonlinear equations. This Special Issue aims at collecting papers presenting new developments in the field of flow and transport in porous media. The 25 published papers deal with different aspects of physical processes and applications such as unsaturated and saturated flow, flow in fractured porous media, landslide, reactive transport, seawater intrusion, and transport within hyporheic zones. Based on their objectives, we classified these papers into four categories: (i) improved numerical methods for flow and mass transport simulation, (ii) looking for reliable models and parameters, (iii) laboratory scale experiments and simulations, and (iv) modeling and simulations for improved process understanding. Current trends on modeling fluid flow and transport processes in porous media are discussed in the conclusion.

1. Introduction

Modeling flow and transport in porous media (FTPM) has become essential in several scientific fields such as chemistry, biology, agriculture, and mechanics, as well as in several engineering disciplines such as petroleum, civil, thermal, and environmental engineering. In water resources, until early 1950, the studies of FTPM were limited to the upper thin layer of soil (few meters) and mainly concerned agriculture applications. Currently, these studies are encountered in a broad spectrum of applications. This covers, among others, development of strategies for groundwater resources management and monitoring [1,2], soil and aquifers contamination [3], predictive studies considering climate change and anthropogenic stresses [4,5], groundwater protection [6], seawater intrusion in coastal aquifers [7], thermal storage in aquifers [8,9,10], management of the vadose zone [11,12], geological sequestration of CO₂ in aquifers [13,14], environmental management of mines [15], surface–subsurface water exchange [16], waste disposal, and geological repository of nuclear waste [17]. In these applications, models of FTPM are mainly used for understanding physical processes and verification of experimental studies. Combined with field observations and uncertainty analysis, they are also used in field studies for practical purposes such as system designing, predictions studies, and making decisions.

However, with these innovative applications, the modeling of FTPM faces specific challenges such as complex physical, chemical, and biological processes occurring at different scales, large length-scale (measured in tens of kilometers or even hundreds of kilometers), large time-scale (measured in tens of years), high heterogeneity, presence of fractures, necessity of 3D simulations, anisotropy, and complex velocity fields. In addition, in real applications, a large number of simulations are often required for model calibration and uncertainty analysis. Model calibration is a crucial step that aims at improving the models realism by confronting simulations to field observations [18]. A reliable modeling analysis requires an uncertainty analysis that aims at understanding how uncertainties in models inputs can propagate to impact the model outputs [19].

Thus, the research methodology in modeling FTPM requires multidisciplinary knowledge and expertise from mathematicians, physicists, chemists, computational scientists, hydrogeologists, and geoscientists. Advances in this topic combine efforts on:

-

Physical and theoretical analysis to improve our understanding of physical processes and key factors that control these processes.

-

Mathematical analysis to establish the equations that can describe these processes and to investigate the solution characteristics of these equations (i.e., solution singularity and multiple solution characteristics).

-

Computational analysis to develop advanced numerical methods, techniques, procedures, and algorithms to improve the accuracy and the efficiency of the numerical models on the basis of contemporary computational science knowledge and expertise.

-

Experimental studies (either at laboratory or field scale) to provide data and information required for process understanding, field characterization, or for the validation of the numerical results.

-

Development of synthetic benchmarks with exact solution (analytical, semianalytical or reference solutions) in order to verify the accuracy and the robustness of the developed numerical models through the comparison between numerical and benchmark solutions.

-

Inverse modeling to estimate indirectly, using appropriate optimization algorithms, model parameters if they are missing or incomplete.

-

Uncertainty and sensitivity analysis in order to estimate the effect of the physical parameters on the model outputs.

In this context, this Special Issue was suggested to collect recent developments in the field of modeling FMTP. A total of 25 papers were published in this issue. The papers address different physical aspects as surface–subsurface flow interaction, seawater intrusion in coastal aquifers, aquifers contamination, flow in unconfined or confined aquifers, reactive transport, flow in fractured porous media, unsaturated flow, landslide, transport within hyporheic zones, and soil internal erosion. For our synthetic analysis in this editorial, we classified the published manuscripts into four categories: (1) improved numerical methods for flow and mass transport simulation (6 papers), (2) looking for reliable models and parameters (8 papers), (3) laboratory-scale experiments and simulations (4 papers), and (4) modeling and simulations for improved process understanding (7 papers). Based on this classification, in the following sections, we present a synthetic analysis of the 25 published manuscripts. A general conclusion and new trends on modeling FTPM are discussed in the last section of this manuscript.

2. Improved Numerical Methods for Flow and Mass Transport Simulation

In the last decade, significant advancements have been achieved in the development of sophisticated numerical models for FTMP. However, due to the broad range of applications, several questions about the reliability, accuracy, and robustness of the numerical models remain open [20]. In addition, there is an ever-growing need for new numerical methods, techniques, and algorithms to take advantage of the new computer designs and continued improvements in computing technology [21].

In this context, Ku et al. [22] considered the solution of the free surface flow problem in layered heterogeneous soils. The problem is treated as a moving boundary problem with a nonlinear boundary conditions. The authors proposed a new method combining an iterative scheme with the collocation Trefftz meshless method. The new method has been verified through several examples and the results show the high accuracy of the proposed method, especially for heterogeneous problems with large permeability contrasts.

In Ku et al. [23], the authors extended the previous method to the solution of subsurface flow problems with a transient moving boundary. The proposed solution is based on the spacetime collocation method with complete Trefftz basis functions. The results obtained show that the proposed method may acquire field solutions with high accuracy and is advantageous for solving inverse problems.

Aharmouch et al. [24] proposed a fully implicit finite volume scheme for efficiently simulating seawater intrusion in coastal aquifers. The proposed model solves the sharp interface nonlinear freshwater–seawater model using Newton’s method. The simulation of a real field test case dealing with the contamination of the Souss–Chtouka aquifer under a large stress period along with high pumping withdrawals shows the robustness of the proposed model.

An integral transformation method to solve the linearized Boussinesq equation that governs the transient groundwater level in a sloping unconfined aquifer is described in Wu and Hsieh [25]. The proposed method is shown to be efficient. Further, it can be directly applied to groundwater level responses under time-varying recharge.

Bahar and Gurarslan [26] addressed a new method of line approach for the solution of the bidimensional advection–dispersion transport equation. The new scheme uses quantic B-spline functions for the spatial discretization with an adaptive Runge–Kutta formula for the time integration. The proposed method is shown to be fast and could be a good alternative for simulation of contaminant transport problems.

Amir and Kern [27] proposed a new method for an efficient solution of the large nonlinear system governing reactive transport in heterogeneous porous media. The system combines partial differential equations describing the transport of chemical species with nonlinear algebraic or differential equations describing the chemical reactions. The authors used a new approach based on the Newton–Krylov method applied to a reduced system of equations. The method was extended to two-dimensional heterogeneous geometries and to mineral precipitation and dissolution reactions. The performances of the method have been shown on synthetic examples.

3. Looking for Reliable Models and Parameters

Fallico et al. [28] conducted an experimental study to investigate the variation modalities of the scaling laws of hydraulic conductivity and effective porosity, according to the variation of the porous medium. A series of slug tests were carried out in a 3D sand box confined aquifer. New relationships, valid for coarse-grained porous media, were determined.

Guérillot and Bruyelle [29] proposed a new approach to upscale transmissibility values. The aim of this procedure is to calculate average reservoir properties from parameters that are highly variable in space (high-resolution geological cells). The proposed method uses the half-block approach and combines the finite volume principles with algebraic methods to provide an upper and a lower bound of the upscaled transmissibility values. Results show that the new approach improves the accuracy of field-scale flow simulation for highly heterogeneous reservoirs and is more accurate and more efficient than the classical transmissibility upscaling methods.

Kozlowski and Ludynia [30] investigated empirical equations to estimate the value of the permeability coefficient for poorly permeable soils. The measurements are carried out on samples of eight poorly permeable soils with different mineralogical composition, origin, and physical and structural properties. A statistical analysis was carried out to identify correlation between soil properties and permeability coefficient. Results show that the three single-variable models that fit the experimental data best are based on the plasticity index, the average pore diameter, and the convexity of silt fraction particles.

Aristodemo et al. [31] performed a series of slug tests in 3D laboratory physical models of different porous media. The authors investigated different filtering methods for recorded hydraulic head data. Fourier and wavelet transforms were applied to eliminate the unwanted high-frequency disturbances in the time variation of the head. A careful smoothing analysis is shown to be very useful when the data sets obtained by slug tests are very large. Further, the mesh wavelet is superior to the other methods (FT and Morl wavelet) for the analyzed statistical parameters.

Yin et al. [32] investigated numerically the fluid flow through rough fractures using the lattice Boltzmann method. The latter was shown to be an effective tool to solve the fluid flow problem with a complex boundary condition at desirable accuracy. Results of simulations showed that the joint roughness coefficient is not an ideal choice to characterize the fluid flow in fractures. An equation with a simple form has been proposed to estimate the permeability based on the aperture and the root mean square of the first derivative of profile of the fracture. The proposed equation has been validated for a synthetic fracture of a wide range of fractal dimensions and standard deviations.

Huang et al. [33] conducted a series of laboratory tests to investigate the permeability of a transparent soil. For such soils, the solid material is replaced by an artificial transparent particle, and the pore fluid by a synthetic liquid with the same refractive index as the particle. The authors introduced a hierarchical approach to quantitatively evaluate the effect of the particle shape on the permeability. Based on the experimental results, a modified Kozeny–Carman model was proposed to estimate the permeability of transparent soils.

Ercan and Kavvas [34] solved transient water flow in porous media in fractional time and multidimensional fractional soil space in anisotropic media by coupling the fractional continuity and motion equations with constitutive relationships. The latter are obtained by least square fits through hydraulic measurements. The results show that fractional governing equations can better mimic the physical process than the traditional Richards equation, as they can capture sub-, super-, and normal-diffusive soil water flow processes during infiltration.

Baalousha et al. [35] investigated the effect of pilot-points number and locations on the calibrated parameters using the pilot-points method. A 3D groundwater flow model was built for the northern karst aquifer of Qatar and calibrated using the parameter estimation and uncertainty analysis (PEST) package. The results revealed that the pilot-points number, locations, and configurations have a significant effect on the calibrated parameter, especially in the high permeable regions corresponding to the karstic zones.

4. Laboratory Scale Experiments and Simulations

Belfort et al. [36] conducted a drainage/imbibition and a rainfall unsaturated flow experiments on a large laboratory experimental tank. A new noninvasive method for water content measurements was proposed. The method, based on photographs, produces 2D water content maps after classical image analysis steps (normalization, filtering, background subtraction, scaling, and calibration). The method does not require any tracer or dye to be injected into the flow tank. The accuracy of the method was demonstrated by comparison against results of numerical simulations performed with an advanced numerical code for solving Richards’ equation.

Isch et al. [37] compared soil hydraulic and transport properties estimated from three bare field plots and in six bare soil lysimeters measurements. A faster transport of bromide was observed in the lysimeters compared to the field plots. Simulation of the bromide transport in lysimeters required the use of the mobile/immobile water model. Results of the optimization procedure showed that most of the estimated parameters for field plot and lysimeter experiments were significantly different. Therefore, parameters optimized for the lysimeters cannot be transposed to field plots.

Kanzari et al. [38] conducted an experiment on a tomato crop cultivated on a silty-clay soil in the semiarid region of Tunisia and irrigated with saline waters of different concentrations. Observed soil salt profiles indicated progressive soil salinization when the plot was irrigated with the most saline water. Observed data were then used to calibrate the Hydrus-1D model for water flow and salt transfer in soils. The model satisfactorily simulated measured yields and reproduced the effects of irrigation with saline waters on crop yields. The model also allows study of the combined effects of soil osmotic and soil matrix potentials on root water uptake.

Younes et al. [39] developed advanced experimental and numerical techniques for an accurate estimation of the soil parameters. Bromide, used as a nonreactive contaminant, was injected at the surface of a large undisturbed soil column and followed by a leaching phase. All hydraulic and mass transport parameters were inferred using a one-step Bayesian estimation with the Markov-chain Monte Carlo sampler. The results prove that the pressure and concentration measurements are able to identify almost all hydraulic and mass transport parameters. The mobile–immobile transport model better reproduces the infiltration experiment. It produces narrower uncertainty intervals for soil parameters and predictive output concentrations.

5. Modeling and Simulations for Improved Process Understanding

Emadi-Tafti and Ataie-Ashtiani [40] investigated causes of landslide triggering. They developed an integrated two-dimensional slope stability model that considers various aspects of hydrological effects and vegetation impacts on the stability of slopes. The analyses of benchmark problems demonstrate that the consideration of matric suction in the unsaturated zone can increase the safety factor more than 90%. It is also observed that the existence of trees with high density on a slope can increase the factor of safety about 50% and prevent shallow landslides.

Liu et al. [41] studied the effect of sediment heterogeneity on conservative and sorptive solute transport within hyporheic zones. Numerical simulations showed that sediment heterogeneity significantly enhances hyporheic exchange and skews solute breakthrough behavior. Strongly heterogeneous sediments with high sorptive capacity compress the solute mixing zone and return more solute mass to the river than homogeneous sediments, even when the bulk physical and chemical properties are equivalent. Permeable, hydraulically connected facies form preferential flow pathways that increase hyporheic exchange and, therefore, solute delivery, but shortens solute breakthrough curves along shallow flow paths.

Lu et al. [42] studied a temporary soil leaching system in which bundled maize straw (straw drainage module) was operated as a subsurface drainage tube and diluted seawater was used for leaching. A numerical model (HYDRUS-3D) based on field measured data was used to simulate the entire leaching process. Results showed that the mass of salt removal was about 1.7 times that of the salt input from the diluted seawater. The salt removal efficiency and soil desalination rate both were negatively related to the seawater mixture rate but were positively associated with the amount of leaching water.

Liu et al. [43] studied dewatering of deep foundation pits excavated in highly permeable media using a finite difference method and showed evident three-dimensional (3D) seepage characteristics. A dewatering design method based on the prediction formulas was proposed in sand and gravel strata and applied to a field case for verification. Measured results show that dewatering is enhanced by 3D flow, forming appropriate pressure distributions for dewatering construction.

Steding et al. [44] investigated leaching zones within potash seams using reactive transport simulations. The authors developed a new reactive transport model complemented by an innovative approach to calculate the interchange of minerals and solution at the water–rock interface. The results of a scenario analysis based on a carnallite-bearing potash seam show that the evolution of leaching zones depends on the mineral composition and dissolution rate of the original salt rock, and that the formation can be classified according to the dimensionless Péclet and Damköhler numbers.

Wang et al. [45] reviewed the mechanisms of soil internal erosion studied by numerical simulation using the discrete element method. The latter, combined with hydromechanical models allowed for modeling the soil deformation characteristic and particle migration, and the evolution of the pore channels. The stability of the soil is mainly determined by force chains and pore channels, whereas the hydraulic conditions act as external disturbances. The erosion process is accompanied by contact failure, force chain bending, kinetic energy burst of particles, and other processes due to multifactor coupling.

Ding et al. [46] revisited the application of fractional flow theory to chemical-enhanced oil recovery (CEOR). The paper provides a review of the analytical modeling of surfactant flooding, polymer flooding, and surfactant–polymer flooding. The effects of various influencing factors on the transport of chemical components and consequently on final incremental oil recovery are investigated, which provide a sharpened insight into designing a successful CEOR process by using the analytical approach.

6. Conclusions

Papers published in this Special Issue provide new theoretical developments for the improvement of model reliability dealing with flow and transport in porous media. This covers development of new numerical techniques that improve accuracy and efficiency. It also covers the development of several analytical or semianalytical solutions that can be helpful in understanding physics and for the assessment of numerical solutions. Beside theoretical developments, the several published papers provide new physical insights by using numerical simulations or by confronting simulations with laboratory experiments. Some papers present new approaches for estimation of model parameters. The diversity of the physical aspects discussed in these papers and the wide range of applications show the significance of flow and transport processes in porous media. These papers make a valuable contribution in improving models and in providing physical insights, but they also pave the way for further research and developments.

Current trends on FTPM concern the improvement of the capacity of the models in representing multiphysical/multiscale processes, and the development of new algorithms and techniques to improve model performances at large scales. For instance, robust simulations of flow and transport processes in fractured porous media at regional scale are currently a real challenge for numerical models. In addition, interaction between different models and integration of new data to improve model realism is also a challenging task. Further, the development of accurate and efficient surrogate FTMP models for sensitivity and uncertainty analyses that require several model evaluations is also a hot topic. Likewise, integration of advanced techniques of machine learning to assist physical models constitutes an important challenge for modeling FTPM.