Acta Polytechnica https://doi.org/10.14311/AP.2025.65.0341 Acta Polytechnica 65(3):341–348, 2025 © 2025 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague SOME REMARKS ON A MATHEMATICAL MODEL FOR WATER FLOW IN POROUS MEDIA WITH COMPETITION BETWEEN TRANSPORT AND DIFFUSION Judita Runcziková, Jan Chleboun, Chiara Gavioli∗, Pavel Krejčí Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mathematics, Thákurova 7, 166 29 Prague, Czech Republic ∗ corresponding author: chiara.gavioli@cvut.cz Abstract. The contribution deals with the mathematical modelling of fluid flow in porous media, in particular water flow in soils, with the aim of describing the competition between transport and diffusion. The analysis is based on a mathematical model developed by B. Detmann, C. Gavioli, and P. Krejčí, in which the effects of gravity are included in a novel way. The model consists of a nonlinear partial differential equation describing both the diffusion and the gravitational transport of water. Although analytical solutions can be obtained for some special cases, only numerical solutions are available in more general situations. The solving algorithm is based on a time discretisation and the finite element method and is written in Matlab. The results of the numerical simulations are shown and the behaviour of the model is discussed. Keywords: Porous media, water flow, transport, diffusion, PDE, numerical methods. 1. Introduction Mathematical models for water flow in unsaturated porous media, such as soils, play an important role in many applications. They have proven to be useful in geoscience and environmental engineering (e.g. spread of contamination into soil and groundwater, and their sanitation), in civil engineering (e.g. the effect of wa- ter content on structures), and in other industrial applications. A simple mathematical model for water flow in saturated porous media is obtained from the mass balance equation (i.e. the continuity equation) com- bined with the Darcy law for the water mass flux. If the porous medium is unsaturated, the water flux is described by a saturation-dependent variant of the Darcy law (also called the Darcy-Buckingham law), and the resulting equation is known as the Richards equation [1]. Many numerical methods have been used to solve the Richards equation. Some of them, with their advantages and limitations, are studied for example in [2–6]. In [7], a comparative analysis is made between vertical water flow models derived from numerical solutions of the Richards equation and reservoir cascade schemes used in hydrological modelling. Additionally, dual-porosity models offer an alternative framework for simulating water flow in unsaturated soils, as exemplified by [8]. While this paper focuses on a macroscopic-scale model examining the competition between transport and diffusion, it is worth noting that the competitive roles of driving forces in porous media have also been explored at the microscopic level, for instance for the pore-network model (discrete model) in [9]. More elaborate mathematical models incorporating additional effects can be found in the literature. Prefe- rential flow (i.e. preferential paths through which the water flows more easily [10]), observed for example in [11], has been addressed in [12] and in the models summarised in [13]. Furthermore, in [14], the Richards equation with hysteresis is introduced. The effects of capillary hysteresis on finger flow, particularly within non-equilibrium Richards equation models, are inves- tigated in [15]. Finally, a semi-continuum model cap- turing both diffusion-like and finger-like flow patterns is shown in [16]. This paper is structured as follows. In Section 2, the partial differential equation (PDE) that models unsaturated fluid flow in a porous medium is presented, and the solving algorithm is described. In Section 3, the results of numerical simulations are shown and discussed. Examples with different initial and boun- dary conditions related to elementary hydrological processes in a soil column are considered in this paper. 2. Materials and methods 2.1. The mathematical model A mathematical model has been proposed by Bettina Detmann, Chiara Gavioli, and Pavel Krejčí in [17] to describe a fluid flow in an unsaturated porous medium under the influence of both diffusion and gravitational transport. Its main feature is a “stickiness” condition: if the saturation stays below a certain threshold s̄, no transport takes place, and water flows only by diffusion. This behaviour is encoded in (1) below, where the notation ( · )+ represents the nonnegative part of a function. The model is a good approximation of the real behaviour when s stays away from the (expected) 341 https://doi.org/10.14311/AP.2025.65.0341 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en J. Runcziková, J. Chleboun, C. Gavioli, P. Krejčí Acta Polytechnica maximum value 1. The full saturation case could be taken into account by including a term of the form (1−s)+ in the transport term. This will be the subject of future research. The mathematical model derived from the mass balance equation is the following: st − κ∆p − 2αg (s − s̄)+ sz = 0, (1) where s = s(x, t) ∈ (0, 1) is the saturation (i.e. the ratio between the water volume and the pore volume), st its time derivative, sz its derivative with respect to z (i.e. the variable representing the vertical direc- tion), and p is the capillary pressure. The time is denoted by t ≥ 0, the porous body is represented by a domain Ω ⊂ R3 with coordinates x = (x, y, z) ∈ Ω and with Lipschitz boundary ∂Ω. The parameter κ is the coefficient, which corresponds to the permeability of the medium, α is the characteristic time related to the friction at the liquid-solid interface, g is the gravity constant, and s̄ ∈ (0, 1) denotes the critical saturation value. If s ≤ s̄, then s can increase or decrease only due to the flow driven by diffusion. The model is complemented by a constitutive relation be- tween the saturation s and the capillary pressure p. Experimental evidence (see e.g. [18] for tests in a close- to-natural environment) indicates that the dependence is of hysteresis type. However, in the case of a purely hysteretic pressure-saturation relationship, the degen- eracy of the problem could affect the convergence of the numerical scheme, see [19]. For simplicity, and as a first step towards more complex models, only linear dependence between s and p is considered in this paper. The model in [17] is also coupled with an initial condition: s(x, 0) = s0(x), (2) and boundary conditions assuming no flow through the vertical parts of the boundary: Q · n = 0, (3) where Q is the water mass flux vector, n the unit outward normal vector, and · denotes the scalar pro- duct. The model is completed by boundary conditions prescribing the inflow/outflow through the horizontal parts of the boundary: Q · n = β(s − sout), (4) where β depends on the permeability coefficient and sout is the saturation outside the domain. In what follows, other boundary conditions on the horizontal part of the boundary will also be taken into account for this model, including the Dirichlet type for its numerical simplicity. The model equation, the initial condition, and the boundary conditions for the 1-D flow are listed in Section 2.2. More details can be found in [17]. 2.2. Algorithm The vertical flow of water through a column of soil is modelled as 1-D flow through an unsaturated porous body. Transport and diffusion effects are taken into account. The relation between pressure and saturation is chosen to be linear, that is: s = γp, (5) where γ = 1. Equation (1) in the 1-D case coupled with Equa- tion (5) reduces to: st − κszz − 2αg (s − s̄)+ sz = 0, (6) where κ > 0, α > 0, g > 0, and s̄ ∈ (0, 1) are constants. The domain is defined as Ω = (−h, 0) with coordinate z ∈ (−h, 0), so that s = s(z, t) and p = p(z, t). The equation-solving algorithm, based on the finite element method, is developed in Matlab. First, pa- rameter values, time interval, and spatial mesh are set up, as well as the initial and boundary conditions. Then, the pdepe procedure (namely, the solver for systems of parabolic and elliptic PDEs of one spatial variable x and time t) is used. The options of the solver are set to their default values, except the relative error tolerance (’RelTol’), which is set equal to 10−5 in the examples in Section 3. The domain is discretised by fixed division points, the spatial mesh is prescribed by matrices of x-coordinates and of y-coordinates and of the numbering of these points. The time mesh is set up automatically, because the time is discretised by the ode15s solver. This multistep solver is a variable- step, variable-order solver based on the numerical differentiation formulas of orders 1 to 5 [20]. More details about the pdepe solver can be found in [21]. Then, the graphical outputs are generated. The initial condition (i.e. the saturation at the ini- tial time t = 0) is defined for all z ∈ (−h, 0) as: s(z, 0) = s0(z), (7) where s0 is a given function in L2(−h, 0), taking values between 0 and 1. The boundary conditions (BC) can be defined for z = −h, 0 and t > 0 as: • given values on the boundary (Dirichlet BC): s(0, t) = f+(t), s(−h, t) = f−(t); (8) • given mass flux through the boundary (Neumann BC): κsz + αg [(s − s̄)+]2 = { f+ for z = 0, f− for z = −h; (9) • mass flux proportional to inner/outer saturation difference (Newton BC): κsz + αg [(s − s̄)+]2 = { −β+(s − s+) for z = 0, β−(s − s−) for z = −h; (10) 342 vol. 65 no. 3/2025 Some remarks on a mathematical model for water flow in porous media where f+, f− are given time dependent functions, β+ > 0, β− > 0 are the permeability coefficients of the top and of the bottom, respectively, and s+, s− ∈ [0, 1] are the saturation values outside the domain in the pos- itive and negative axis directions, respectively. BC, as in (8) with f+ = f− = 0, is considered in Section 3.3, BC, as in (9) with f+ = f− = 0, is considered in Section 3.1 and Section 3.2. The validity of the mass conservation is monitored by computing the definite integral of the solution s over the domain by trapezoidal numerical integration. In the case of no inflow/outflow, the value of the integral has to remain constant over time. In the other cases, the value of the integral should comply with the difference between the inflow and outflow of the water through the boundary. Thus, the water mass balance can be checked during the execution of the Matlab program. 3. Results An example of redistribution after infiltration is pre- sented in Section 3.1. The process of redistribution is described in [22] and is also mentioned in [23], where hysteresis effects are additionally included. The ini- tial condition corresponds to saturated soil in the top part of the column, within a given depth below the terrain, with a steep transition (almost a jump) to dry soil. In Section 3.2, the process of wetting the soil column with a water source located at the bottom of the column is examined. In Section 3.3, the effect of κ and s̄ on a similar setting is studied. 3.1. Example 1 In this first example, representing the redistribution after infiltration, the parameters are chosen as follows: permeability κ = 0.005, 2αg = 1, s̄ = 0.25, maximal depth h = 5, value of the space discretisation parame- ter d = 0.01. The value of the time step τ is chosen automatically by the solver. The initial condition (displayed in Figure 1) is chosen as: s0(z) =  1 for z > −0.50, 100z + 51 for z = [−0.51, −0.50], 0 for z < −0.51, (11) which means saturated soil to the depth of 0.5 and dry soil below. The BC is that of impermeable boundary: κsz + αg [ (s − s̄)+]2 = 0, (12) so that the control integral of the saturation over the domain remains constant over time, as displayed in Figure 2. The only exception is the small interval (of order 10−4) at the beginning, caused by a sharp change in the shape of the saturation curve. The saturation profile at selected times t = 0.5, 5, 250 and 2 500 is plotted in Figure 3. At the beginning of the redistribution process, a rapid change in the saturation profile in the upper part can be noticed. Figure 1. Initial condition (11). Figure 2. Integral of the solution over the domain. A significant shift of the water content from the upper to the lower part is then observed, that is, the water content decreases in the upper part (above the wetting front) and increases in the lower part. The terrain saturation decreases from s = 1 to near s = s̄. This is caused by the dominant effect of transport in this zone. The effect of diffusion (note that κ has a significantly large value) can be seen on the shape of the wetting front, which is slowly tilting and moving downwards. The saturation value over time decreases close to the terrain and increases in the lower part. Accord- ing to (6), if the value s(zi, t) of the saturation at a given depth zi is smaller than s̄, then the transport term is not active at zi, and the water movement is driven only by diffusion. Hence, if the saturation in the upper part decreases until it reaches the value s̄, then only the diffusion term contributes to the redis- tribution of water. With the above parameters, the saturation along the entire soil column drops below s̄ at time t = 129, as can be seen in Figure 4. The water is then redistributed by diffusion until it reaches a uniform distribution profile along the soil column, and the maximal value of s decreases until it meets the minimal value, as shown in Figure 4. From time 343 J. Runcziková, J. Chleboun, C. Gavioli, P. Krejčí Acta Polytechnica Figure 3. Evolution of the solution s(z, t) over time with κ = 0.005 and initial condition as in (11). Figure 4. Minimum and maximum of the solution s(z, t) in the range defined by the time steps and the critical saturation value. t = 609, the difference between the maximum and minimum value of the solution is less than 0.1. The process slows down as the saturation value decreases. 3.2. Example 2 The initial condition is now, i.e. in the example of soil wetting from the bottom water source, changed to: s0(z) =  0 for z > −4.50, −30z − 135 for z = [−4.50, −4.51], 0.3 for z < −4.51, (13) which means unsaturated zone with a thickness of 0.49 with saturation s = 0.3 (hence s > s̄) in the bottom part, with a steep transition to dry soil up to the terrain. The other parameters and boundary con- ditions are the same as above. The saturation profile at selected times t = 0.5, 5, 250 and 2 500 is plotted in Figure 5. First, the flow driven by transport near the bottom and by diffusion above the unsaturated part is observed. Then the transport-driven flow in the lower part becomes weaker, the saturation starts to decrease, and when the threshold value s̄ is met, the transport term becomes inactive, so that only the upper part are continuously wetted by diffusion. The profiles of minimal and maximal soil saturation over time are shown in Figure 6, where the peak in the bottom figure occurs at the time when the increase in saturation due to transport in the lower part ends. 344 vol. 65 no. 3/2025 Some remarks on a mathematical model for water flow in porous media Figure 5. Evolution of the solution s(z, t) over time with κ = 0.005 and initial condition as in (13). 3.3. Example 3 The influence of the choice of κ and s̄ on the saturation profile is shown in the next examples. The following parameters are chosen: γ = 1, 2αg = 1, permeability κ in the range [0, 0.01], critical saturation value s̄ = 0, s̄ = 0.25 and s̄ = 0.40, maximal depth h = 5, and the value of the space discretisation parameter d = 0.01. The value of the time step τ is chosen automatically by the solver. The initial condition is now chosen as s0(z) =  −0.5z for z > −2.00, 100z + 201 for z = [−2.00, −2.01], 0 for z < −2.01, (14) and is displayed in Figure 7. For numerical simplicity, the boundary condition here is of Dirichlet type (8) with f+ = f− = 0. The diffusion-driven flow is determined by the value of the permeability κ. If κ = 0, the diffusion term in (6) is inactive. If the value of κ is close to 0, for example κ = 0.001, the diffusion is negligible. Higher values of κ are associated with strong diffusion. If the diffusion is strong enough, that is, the value of κ is large enough, it is possible to choose an initial condition with a steep transition between two different saturation values, because the gradient of the solution soon becomes slightly smaller during the calculation due to the diffusion flow, and the simulation is stable. This situation occurs in this example in Figure 8, bottom, and also in Sections 3.1 and 3.2. Conversely, when the value of κ is small, the gradient of the solution remains large, or even tends to infinity, and numerical instabilities (i.e. oscillations of the solution) occur. Depending on the value of κ, these oscillations can cause calculation collapse, as shown in Figure 8, top, or disappear, as in Figure 8, centre. The same analysis applies if a high saturation gra- dient between the wet and the dry parts of the porous body occurs later in the calculation. In Figure 9 the effect of the critical saturation value s̄ on the solution s is shown. When s̄ is higher, the wetting front reaches a smaller depth. 4. Conclusions A mathematical model for the unsaturated flow of water in a porous medium was studied. Numerical simulations of a vertical 1-D water flow in an unsatu- rated soil column were performed, and the effect of two parameters, the permeability κ and the critical saturation value s̄, was studied. Their magnitude has a significant impact on the solution and, in some cases, on the stability of the calculation. 345 J. Runcziková, J. Chleboun, C. Gavioli, P. Krejčí Acta Polytechnica Figure 6. Minimum and maximum of the solution s(z, t) in the range defined by the time steps (top) and zoom to the beginning of the calculation (bottom). Figure 7. Initial condition (14) with a steep wetting front. The case of redistribution of water in soil after in- filtration was shown and the results were compared with theoretical physical expectations from the tech- nical literature. The resulting saturation profiles are in agreement with [22]. The case of wetting from below due to pressure difference was also considered. The model has shown a good approximation of the physical behaviour because, with the no-flux boun- dary condition (12), the saturation values stay under control even if we initially allow s = 1 in some zone of given thickness as in (11). Initially, the model can be used for simulating lab- oratory experiments, such as water flow through soil (a). Without diffusion (κ = 0). (b). Negligible diffusion (κ = 0.001). (c). Strong diffusion (κ = 0.01). Figure 8. Effect of permeability κ on the solution s(z, 0.5), with s̄ = 0.25. columns. For application to regional water flow models within catchment areas, more complex model formula- tions are required. Furthermore, the influence of hys- teresis, which can be incorporated using the Preisach model to describe the relationship between satura- 346 vol. 65 no. 3/2025 Some remarks on a mathematical model for water flow in porous media (a). Without stickiness (s̄ = 0). (b). With lower stickiness (s̄ = 0.25). (c). With higher stickiness (s̄ = 0.4). Figure 9. Effect of s̄ on the solution s(z, 5) with κ = 0.005. tion and capillary pressure (similar to the approach in [19]), and preferential flow should be thoroughly investigated. From the modelling point of view, future research will include the case where the saturation s ap- proaches 1. From the numerical point of view, the investigation will be directed towards the development of an algorithm for a 2-D flow, to test the behaviour of the model in 2-D (Ω ⊂ R2). List of symbols s Saturation p Capillary pressure κ Permeability s̄ Critical saturation value α Characteristic time g Gravity constant Acknowledgements Supported by the project Centre of Advanced Applied Sciences (CAAS) No. CZ.02.1.01/0.0/0.0/16 019/0000778, which is co-financed by the European Union; by the Grant Agency of the Czech Technical University in Prague, grant No. SGS23/092/OHK1/2T/11; by the European Union’s Horizon Europe research and innovation pro- gramme under the Marie Skłodowska-Curie grant agree- ment No. 101102708; by the Czech Science Foundation (GAČR) project No. 24-10586S. References [1] R. Černý, P. Rovnaníková. Transport processes in concrete. Spon Press, London, 2002. ISBN 0-415-24264-9. [2] R. Haverkamp, M. Vauclin. A comparative study of three forms of the Richard equation used for predicting one-dimensional infiltration in unsaturated soil. Soil Science Society of America Journal 45(1):13–20, 1981. https://doi.org/10.2136/sssaj1981. 03615995004500010003x [3] P. J. Ross. Efficient numerical methods for infiltration using Richards’ equation. Water Resources Research 26(2):279–290, 1990. https://doi.org/10.1029/WR026i002p00279 [4] F. List, F. A. Radu. A study on iterative methods for solving Richards’ equation. Computational Geosciences 20(2):341–353, 2016. https://doi.org/10.1007/s10596-016-9566-3 [5] M. W. Farthing, F. L. Ogden. Numerical solution of Richards’ equation: A review of advances and challenges. Soil Science Society of America Journal 81(6):1257–1269, 2017. https://doi.org/10.2136/sssaj2017.02.0058 [6] H. Kamil, A. Beljadid, A. Soulaïmani, Y. Bourgault. Semi-implicit schemes for modeling water flow and solute transport in unsaturated soils. Advances in Water Resources 193:104835, 2024. https://doi.org/10.1016/j.advwatres.2024.104835 [7] C. Gandolfi, A. Facchi, D. Maggi. Comparison of 1D models of water flow in unsaturated soils. Environmental Modelling and Software 21(12):1759–1764, 2006. https://doi.org/10.1016/j.envsoft.2006.04.004 [8] J. Lewandowska, A. Szymkiewicz, K. Burzyński, M. Vauclin. Modeling of unsaturated water flow in double-porosity soils by the homogenization approach. Advances in Water Resources 27(3):283–296, 2004. https://doi.org/10.1016/j.advwatres.2003.12.004 347 https://doi.org/10.2136/sssaj1981.03615995004500010003x https://doi.org/10.2136/sssaj1981.03615995004500010003x https://doi.org/10.1029/WR026i002p00279 https://doi.org/10.1007/s10596-016-9566-3 https://doi.org/10.2136/sssaj2017.02.0058 https://doi.org/10.1016/j.advwatres.2024.104835 https://doi.org/10.1016/j.envsoft.2006.04.004 https://doi.org/10.1016/j.advwatres.2003.12.004 J. Runcziková, J. Chleboun, C. Gavioli, P. Krejčí Acta Polytechnica [9] A. Mohammadi, M. R. Rasaei, V. Mashayekhizadeh, A. Nakhaee. Pore-network modelling of combined molecular diffusion and gravity drainage mechanisms in a porous matrix block: The competitive role of driving forces. The Canadian Journal of Chemical Engineering 101(5):2923–2947, 2023. https://doi.org/10.1002/cjce.24647 [10] S. E. Allaire, S. Roulier, A. J. Cessna. Quantifying preferential flow in soils: A review of different techniques. Journal of Hydrology 378(1–2):179–204, 2009. https://doi.org/10.1016/j.jhydrol.2009.08.013 [11] M. S. Babel, R. Loof, A. Das Gupta. Fingered preferential flow in unsaturated homogeneous coarse sands. Hydrological Sciences Journal 40(1):1–17, 1995. https://doi.org/10.1080/02626669509491387 [12] T. Vogel, J. Březina, M. Dohnal, J. Dušek. Physical and numerical coupling in dual-continuum modeling of preferential flow. Vadose Zone Journal 9(2):260–267, 2010. https://doi.org/10.2136/vzj2009.0091 [13] H. H. Gerke. Preferential flow descriptions for structured soils. Journal of Plant Nutrition and Soil Science 169(3):382–400, 2006. https://doi.org/10.1002/jpln.200521955 [14] B. Schweizer. Instability of gravity wetting fronts for Richards equations with hysteresis. Interfaces and Free Boundaries 14(1):37–64, 2012. https://doi.org/10.4171/IFB/273 [15] W. J. Roche, K. Murphy, D. P. Flynn. Modelling preferential flow through unsaturated porous media with the Preisach model and an extended Richards equation to capture hysteresis and relaxation behaviour. Journal of Physics: Conference Series 1730(1):012002, 2021. https://doi.org/10.1088/1742-6596/1730/1/012002 [16] J. Kmec, M. Šír. Modeling 2D gravity-driven flow in unsaturated porous media for different infiltration rates. Hydrology and Earth System Sciences 28(22):4947–4970, 2024. https://doi.org/10.5194/hess-28-4947-2024 [17] B. Detmann, C. Gavioli, P. Krejčí. Soil moisture transport and diffusion driven by gravity and capillarity. [In preparation]. [18] A. Poulovassilis. Hysteresis of pore water, an application of the concept of independent domains. Soil Science 93(6):405–412, 1962. [19] C. Gavioli, P. Krejčí. Degenerate diffusion with Preisach hysteresis. Discrete and Continuous Dynamical Systems - S 16(12):3677–3708, 2023. https://doi.org/10.3934/dcdss.2023154 [20] L. F. Shampine, M. W. Reichelt. The MATLAB ODE suite. SIAM Journal on Scientific Computing 18(1):1–22, 1997. https://doi.org/10.1137/S1064827594276424 [21] R. D. Skeel, M. Berzins. A method for the spatial discretization of parabolic equations in one space variable. SIAM Journal on Scientific and Statistical Computing 11(1):1–32, 1990. https://doi.org/10.1137/0911001 [22] M. Kutílek, V. Kuráž, M. Císlerová. Hydropedologie 10 [In Czech; Hydropedology 10]. Second printing. CTU in Prague, Prague, 2000. ISBN 80-01-02237-4. [23] E. G. Youngs. Application of scaling to soil-water movement considering hysteresis. In Scaling in Soil Physics: Principles and Applications, chap. 3. John Wiley & Sons, Ltd, 1990. https://doi.org/10.2136/sssaspecpub25.c3 348 https://doi.org/10.1002/cjce.24647 https://doi.org/10.1016/j.jhydrol.2009.08.013 https://doi.org/10.1080/02626669509491387 https://doi.org/10.2136/vzj2009.0091 https://doi.org/10.1002/jpln.200521955 https://doi.org/10.4171/IFB/273 https://doi.org/10.1088/1742-6596/1730/1/012002 https://doi.org/10.5194/hess-28-4947-2024 https://doi.org/10.3934/dcdss.2023154 https://doi.org/10.1137/S1064827594276424 https://doi.org/10.1137/0911001 https://doi.org/10.2136/sssaspecpub25.c3 Acta Polytechnica 65(3):341–348, 2025 1 Introduction 2 Materials and methods 2.1 The mathematical model 2.2 Algorithm 3 Results 3.1 Example 1 3.2 Example 2 3.3 Example 3 4 Conclusions List of symbols Acknowledgements References