Numerical Simulations on the Front Motion of Water Permeation into Anisotropic Porous Media

Water permeation into a porous medium is a common but important phenomenon in many engineering ﬁ elds such as hydraulic fracturing. The water permeation front moves with time and may signi ﬁ cantly impact the ﬁ eld variable evolution near the water front. Many algorithms have been developed to calculate this water front motion, but few numerical algorithms have been available to calculate the water front motion in anisotropic ﬂ uid-solid couplings with high computational e ﬃ ciency. In this study, a numerical model is proposed to investigate the front motion of water permeation into an anisotropic porous medium. This model fully couples the mechanical deformation, ﬂ uid ﬂ ow, and water front motion. The water front motion is calculated based on a directional Darcy ’ s ﬂ ow in the anisotropic porous medium, and a revised formula with a correction coe ﬃ cient is developed for the estimation of permeation depth. After veri ﬁ cation with three sets of experimental data, this model is used to numerically investigate the impacts of permeability, viscosity, permeability anisotropy, and mechanical anisotropy on water front motion. Numerical results show that the proposed model can well describe the anisotropic water permeation process with reasonable accuracy. The permeation depth increases with permeability, mobility, and mechanical anisotropy but decreases with viscosity and permeability anisotropy. The correction coe ﬃ cient mainly depends on porosity evolution, ﬂ ow pattern, mobility, permeability anisotropy, and mechanical anisotropy.


Introduction
Water permeation into a porous medium is a common phenomenon in many engineering fields: masonry structure [1][2][3][4], hydraulic fracturing [5][6][7][8][9], concrete [10][11][12][13], and carbon capture and storage (CCS) [14][15][16]. This phenomenon may cause the frost and salt damage to brickwork in masonry structures and the rot damage to wood structures [1]. It may also cause the corrosion of building surface by the wind-driven rain [2] and thus significantly affect the durability of masonry structures [3]. In hydraulic fracturing, water permeation has significant effects on the stress redistribution and fracture behaviors of rocks [5,8,17]. During the injection of fracturing fluid, the permeation zone gradually extends with water front motion and the pore pressure increases within water permeation zone. This induces the redistribution of effective stress and affects the initiation and propagation of fractures [18][19][20]. Concrete, as a type of porous medium, is easily penetrated by water. Water has remarkable effects on the properties of concrete. For example, water permeation increases the elastic modulus but decreases the compressive strength and durability of concrete [10,12,21]. The sound velocity and frequency of maximum transmission are affected by the water content of concrete [22]. Besides, in the CCS technology, acidic fluid intrusion has become a potential risk to the safety of geological storage of carbon dioxide (CO 2 ) in underground formations [15,23]. The acidic fluid permeation through the caprock layer affects the groundwater system and reduces the local air quality on the earth surface [24]. From the above, water permeation into a porous medium has become a serious but interesting topic to engineering fields. Therefore, it is important to investigate this phenomenon of water permeation.
During water permeation, water front motion has been focused. Several sets of experiments have been conducted to investigate the permeation depth of water front [4, 11-13, 25, 26]. It is discovered that permeation depth varies with the properties of porous medium, the viscosity of injection fluid, injection time, and injection pressure. Khatri and Sirivivatnanon [25] performed laboratory tests on concrete to determine its permeability based on the permeation depth of water front. Murata et al. [11] experimentally and theoretically studied the evolution of permeation depth with water pressure, pressurized time, and concrete type. Based on water permeation tests, Yoo et al. [12] found that the fluid flow in concrete obeys Darcy's law under low water pressure and the diffusion flow under high water pressure. Al-Maamori et al. [26] conducted fluid permeation tests on Queenston shale by different types of fluid. Their results indicated that permeation depth varies with fluid type. By means of electrical conductivity sensors, Guizzardi et al. [4] investigated the relationship between permeation depth and arrival time in building masonry materials. However, most of these experiments were restricted to samples with simple geometry in 1D case. The anisotropy of porous medium has been ignored, although it is a distinguished feature of porous medium such as rock materials. Thus, it is necessary to extend the study from 1D to 2D and to investigate the effect of anisotropy of porous medium on water front motion.
In numerical simulations, the moving interface, including water front, has crucial effects on fluid dynamics, heat transfer, and other physical processes [27,28]. The boundary shape and field variables significantly change with water front motion. Therefore, the calculation of the permeation depth of water front which represents the moving boundary in water permeation becomes an essential job. This moving boundary problem has been implemented through a fully coupled model which considers the mechanical deformation, the change of pore pressure, and a two-phase flow model [29][30][31]. However, this fully coupled model is so complicated in implementation and so slow in computation. It is also difficult to obtain the real-time permeation depth of water front during the computational process. Several numerical algorithms have been proposed to calculate the permeation depth of water front [13,[32][33][34][35]. Lockington et al. [32] developed a new approach based on Boltzmann transformation to analyze the water permeation profile under exponential diffusivities. Zhou [34] established a general method to solve hydraulic diffusivity from sorptivity tests and determined an approximate solution of Boltzmann variable for any distribution of diffusivity. Li et al. [13] proposed a three-phase mesoscale model to simulate the unsaturated transport in concrete. The permeation depth obtained from the numerical modeling showed a reasonable consistence with experimental results. Further, based on Darcy's law, a simple approach was presented for the estimation of permeation depth with injection time, pressure difference, and pressure magnitude [11,12,15]. A fitting equation was proposed to describe the correlation between permeation depth and injection time for different injection fluids [26]. However, most of these numerical models were established for 1D geometry and ignored the effect of anisotropy on water permeation. The interaction between mechanical behaviors and water permeation was disregarded because the application of these models to the fluid-solid coupling model was so complicated. Moreover, the impacts of permeability, viscosity, permeability anisotropy, and mechanical anisotropy on permeation depth have not been well studied by numerical simulation. Therefore, it is essential to propose a simple and effective method to capture the permeation depth of water front, which should be easily incorporated into the anisotropic fluid-solid coupling model.
In this study, a moving boundary method is proposed to investigate the front motion of anisotropic water permeation. Based on this method, a new anisotropic fully coupled model is established by creatively incorporating water front motion into the interaction between mechanical deformation and fluid flow in an anisotropic porous medium. This model is verified by three sets of experimental data and then used to investigate the impacts of permeability, viscosity, permeability anisotropy, and mechanical anisotropy on the permeation depth of water front. The rest of this paper is organized as follows: Section 2 gives the governing equations for the anisotropic mechanical deformation, anisotropic fluid flow, and constitutive laws for directional permeability as well as the moving boundary method. Section 3 verifies the numerical model with three sets of experimental data. Section 4 discusses the revised formula for the estimation of water front motion. Section 5 conducts a parametric study on permeation depth. The last section summarizes the understandings of this study.

Governing Equations for Each Physical Process
This study made the following assumptions: (1) porous medium is continuous and anisotropic. Its deformation is linearly elastic and infinitesimal; (2) the fluid flow in an anisotropic porous medium follows Darcy's law; (3) the chemical effect of water within porous medium is ignored; (4) all processes involved are isothermal.

Governing Equation for Anisotropic Mechanical
Deformation. For a pseudostatic deformation process, the equation of motion is where σ e ij is the effective stress component; f i is the body force per unit volume in the ith direction; p is the pore pressure; α is the Biot coefficient and α = 1 − K/K s ; K is the bulk modulus of the porous medium; K s is the bulk modulus of grains.

Geofluids
The effective stress is expressed as where σ ij is the total stress component and δ ij is the Kronecker delta.
For an orthotropic linear elastic porous medium, its constitutive model is expressed as where E i is the elastic modulus in the ith direction; ν ij is the Poisson ratio that characterizes the strain response in the ith direction to stress acting in the jth direction; G ij is the shear modulus in the jth direction whose normal is in the ith direction.

Governing Equation for
Anisotropic Fluid Flow. Based on Darcy's law and the mass conservation law, the governing equation for anisotropic fluid flow is expressed as [8] where ϕ is the current porosity; k i is the permeability in the ith direction; S is the current effective volumetric strain; ε v is the current volumetric strain; c ρ is the isothermal coefficient of compressibility which is defined as [36] where ρ is the fluid density. For water, the density is usually constant and c ρ = 0.
Further, a general porosity model for porous medium is obtained as [37] where R m = α/ϕ 0 ; ϕ 0 is the initial porosity; Δε e is the change of effective volumetric strain and expressed as Ignoring the effect of sorption and chemical reaction, the current and initial effective volumetric strains are defined as where ε v0 is the initial volumetric strain; p 0 is the initial pore pressure. Furthermore, the directional permeability in a two-dimensional domain is associated with the effective strain in its perpendicular direction as [8] where k i is the initial permeability in the ith direction; R j = △ε mj /△ε j is the ratio of matrix strain to the strain of whole element in the jth direction.

Moving Boundary Algorithm for Anisotropic Water Front
Motion. The seepage zone expands continuously as the water front moves forward with the injection of fluid [11-13, 26, 38]. This water front motion forms a moving boundary problem in a porous medium. Figure 1 shows a typical 1D seepage-controlled model with moving water front. In this model, the water front moves with the injection of fluid and forms a changing seepage zone. Based on Darcy's law, the moving speed of water front in the ith direction is described as Equation (11) can be used to simulate the anisotropic front motion of water permeation. In this equation, the moving speed is associated with the permeability, viscosity, and pore pressure gradient. The increase of permeability and pore pressure gradient will enhance the water front motion, while the increase of viscosity will hinder it.

Numerical Computation
Procedure. The numerical model proposed in this study is composed of the governing 3 Geofluids equation for anisotropic mechanical deformation (equations (1)-(4)), the governing equation for anisotropic fluid flow (equation (5)), and the moving boundary algorithm for anisotropic water front motion (equation (11)). Their interactions are shown in Figure 2. After incorporating the constitutive laws of equations (7) and (10), this model is implemented by COMSOL Multiphysics for water permeation into an anisotropic porous medium. Particularly, the mechanical deformation equations are solved by Solid Mechanics module, and fluid flow equations are implemented by PDE module. The moving boundary method is conducted by Deformed Geometry module. Through the coupling solver in COMSOL Multiphysics, the proposed numerical model can be solved if the boundary and initial conditions, parameters, and computation time are given.

Verification of Numerical Model
The numerical model is verified by three sets of water permeation tests: two sets are on concrete with different injection pressures and one set is on Queenston shale with different types of injection fluid.
3.1. Verification with Experimental Data by Murata et al. [11]. A set of water permeation tests was conducted to investigate the water tightness of concrete [11]. Three water-cement ratios were set as 0.55, 0.70, and 0.80. Each specimen was 150 mm in diameter and 150 mm in height. The curved surface was maintained watertight, and a specified water pressure was applied to the top end. Photographs were used to measure the permeation depth of water with injection time. In this study, the concrete is simplified as a 2D geometry model with a square of 150 mm × 150 mm. The main computational parameters are listed in Table 1. The injection time is 48 h. Figure 3 compares the numerical results by our model with the experimental data by Murata et al. [11]. This figure shows that the permeation depth by our model agrees well with the experimental data. With the increase of water pressure, the permeation depth increases rapidly at first and its increasing rate decreases gradually. It indicates that the numerical model proposed in this study is effective for the front motion of water permeation into an anisotropic porous medium. Therefore, our model is able to simulate the front motion of 1D water permeation.

Verification with Experimental Data by Yoo et al. [12].
Yoo et al. [12] conducted another set of water permeation tests to investigate the water tightness of concrete. Their water-cement ratios were 0.40, 0.50, and 0.60, respectively. Cylindrical concrete specimens (150 mm in diameter and 300 mm in length) were cast with a hollow cylindrical core of 20 mm in diameter at the center. The water was injected into the inner borehole by controlling time and pressure. In this study, the concrete is simplified as a 2D geometry model of circular ring with 150 mm outer diameter and 20 mm inner diameter. The main computational parameters are listed in Table 2. The water pressure is 1.5 MPa. Figure 4 compares the numerical results by our model with the experimental data by Yoo et al. [12]. From this figure, the permeation depth by our model is in good agreement with the experimental data. With the increase of injection time, the permeation depth increases rapidly at first and its increase rate decreases gradually. Due to the different permeability, the permeation depth varies with the water-cement ratio at the same injection time. Therefore, our model is still capable of calculating the front motion of 2D water permeation.

Verification with Experimental
Data by Al-Maamori et al. [26]. A test procedure was developed to investigate the difference of permeation depth between lubricant fluids and water in Queenston shale. The cylindrical specimens were cut into an approximate 0.24 m in length and 0.063 m in diameter. Fluid pressure was applied on the bottom of the cylindrical specimens with water, bentonite solution, and polymer solution. The pressure increased gradually over 1 h until it reached 200 kPa and then kept constant throughout the test period. In this study, the concrete is simplified as a 2D geometry model with a rectangle of 0 24 m × 0 063 m. The main computational parameters are listed in Table 3. Figure 5 shows the numerical results by our model with the experimental data by Al-Maamori et al. [26]. It shows that the permeation depth by our model well fits the experimental data. The permeation depth increases with injection time. At the same injection time, the permeation depth varies with the type of injection fluid. This is caused by the different viscosity among water, bentonite solution, and polymer solution. Therefore, our model has the capability to determine the front motion of 1D permeation by different types of fluids.

Revised Formula for the Estimation of Water Front Motion
The above section has established a numerical model to investigate the front motion of water permeation into an anisotropic porous medium. This model was well verified with three sets of experimental data. This section will further estimate the water front motion through a revised formula. This is because a simple and effective formula without any

A Revised Analytical Formula for the Prediction of Permeation
Depth. In 1D water permeation, the fluid flow velocity keeps constant over time and space, and its pressure gradient is constant. Ignoring the time-dependent term in equation (5) yields With constant injection pressure and initial pore pressure, the moving speed of water front is expressed as Integrating equation (13) yields The initial condition is L i = 0 when τ = 0. The permeation depth of 1D water permeation is then obtained as [15] where L W is the permeation depth of water front; p 1 is the injection pressure. This formula is based on 1D linear flow and the time-depended term in equation (5) is ignored. It did not consider the anisotropy of porous medium, the fluid-solid coupling effects and other flow patterns. In this study, a    5 Geofluids revised formula is proposed to calculate the permeation depth through introducing a correction coefficient into the 1D formula of equation (15) as where L R is the revised permeation depth; f is a correction coefficient. The correction coefficient in equation (16) comprehensively expresses the impacts of flow pattern, the time-depended term, the anisotropy of porous medium, and the fluid-solid coupling on the permeation depth of water front. Generally, this correction coefficient is so complicated that an accurate theoretical solution is difficult to be obtained if all these influencing factors are considered. It can be determined through fitting numerical solutions or experimental data. For both flow pattern and time-dependent term, following two sections will analytically discuss their impacts on this coefficient.

Flow Pattern Induced Revision of Permeation Depth.
For the radial flow in 2D water permeation, the profile of pore pressure is [39] where R is the radius of the borehole; R L is the radius of permeation depth and R L = R + L i . Substituting equation (17) into equation (11) yields in which With initial condition of L i = 0 when t = 0, the permeation depth can be integrated as Equation (19) shows that the moving speed of water front is smaller in 2D water permeation than in 1D water     6 Geofluids permeation. The correction coefficient f 1 in equation (20) reflects the effect of flow pattern on permeation depth. It may induce obvious errors of permeation depth if the 1D formula is used to predict the front motion of 2D water permeation.

Porosity Evolution Induced Revision of Permeation
Depth. With consideration of time-dependent term, equation (5) can be rewritten as The right term of equation (21) is a function of time and position, which is here defined as Ignoring the evolution of k i can obtain Therefore where C 1 and C 2 are the constants; F x i = βdx i dx i is a function of porosity evolution.
Its boundary conditions are With boundary conditions of equation (25), equation (24) can be solved as Substituting equation (26) to equation (11) yields in which With initial condition of L i = 0 when t = 0, the permeation depth can be integrated as Equation (28) shows that porosity evolution has some impacts on the moving speed of water front. The correction coefficient f 2 in equation (29) reflects the effect of porosity evolution on permeation depth. If the time-dependent item is ignored, the porosity evolution is zero. Thus, F x i = 0 and λ 2 = 1. Equation (29) is equal to equation (15). Otherwise, with the consideration of porosity evolution, F x i ≠ 0. It may decrease the prediction accuracy of permeation depth without the correction coefficient f 2 .
Both equations (20) and (29) show that 1D formula has inevitable limits on calculating the water front motion. Their correction coefficients are so complicated to be theoretically derived. Therefore, the revised formula of equation (16) is essential, and the determination of the correction coefficient will be discussed in the next section.

Parametric Study on Permeation Depth
The above results show that the proposed numerical model in this study can investigate the front motion of water permeation into an anisotropic porous medium. This section will explore key influencing factors on permeation depth through parametric study on permeability, viscosity, permeability anisotropy, and mechanical anisotropy. Their effects on correction coefficients are discussed. Figure 6 describes a geometry model with a square of 1m × 1m and a borehole of 0.1 m in diameter. The injection pressure is taken as 20 MPa, and the initial pore pressure is 5 MPa. The confining pressure is 10 MPa in both the xth and yth directions. The computational parameters are listed in Table 4.

Impact of Mobility
5.1.1. Impact of Permeability. In order to investigate the impact of permeability, the permeability is taken as a range from 0.01μD to 500μD and the viscosity is fixed as 1 mPa·s, 10 mPa·s to 100 mPa·s, respectively. Figure 7 represents the change of permeation depth with permeability at the injection time of 100 h. The permeation depth and its change rate increase with permeability at the same viscosity, while the permeation depth at the same permeability decreases with fluid viscosity. The curve shape of permeation depth versus permeability is similar at different viscosity. This indicates that larger viscosity obstructs water front motion and reduces the permeation depth. Furthermore, the permeation depth by 1D formula and the revised formula are also presented in Figure 7. It shows that the revised formula is more suitable than 1D formula to calculate the front motion of 2D water permeation. The correction coefficient is a function of permeability as Therefore, viscosity only affects the constant term.

Impact of Viscosity.
The viscosity is taken as a range from 0.1 mPa·s to 500 mPa·s when the permeability in the yth direction is fixed as 0.1μD, 1μD, and 10μD, respectively. Figure 8 represents the change of permeation depth with viscosity at the injection time of 100 h. The permeation depth and its change rate decrease with the viscosity at the same permeability, while the larger permeability promotes the water permeation and increases the permeation depth at the same viscosity. The curve shape of permeation depth versus viscosity is also similar among different permeabilities. In addition, the permeation depth by 1D formula is larger than the numerical results. The revised formula can well fit the numerical results. The correction coefficient is a function of viscosity as Obviously, permeability will affect the permeation depth, and this effect varies with fluid viscosity.

Geofluids
where κ i is the directional mobility. Figure 9 shows the change of permeation depth with mobility. All the data are taken from Figures 7 and 8. It is found that the same mobility causes the same permeation depth, even though the viscosity and the permeability are different. The permeation depth increases with mobility, which means that larger mobility may enhance water permeation and the permeation depth is bigger. Also, the permeation depths by 1D formula and the revised formula are also compared in Figure 9. The revised formula is more accurate than 1D formula. The mobility combines the impacts of permeability and viscosity, thus equations (30) and (31) can be unified by a function of mobility as f κ y = −0 072 lg κ y − 0 365 33 5.2. Impact of Permeability Anisotropy. For convenience, the permeability anisotropy is defined as the anisotropic ratio of the permeability in the xth direction (horizontal) to the permeability in the yth direction (vertical): The viscosity is taken as 1 mPa·s, and the permeability in the yth direction is fixed as 5 × 10 -19 m 2 , 1 × 10 -18 m 2 , and 2 × 10 -18 m 2 , respectively. The permeability anisotropy is taken in a range from 0.1 to 10. The permeation depth at the permeability anisotropy of 5 is used as the base, thus the permeation depth ratio is defined as Figure 10 presents the changes of permeation depth and permeation depth ratio with permeability anisotropy at different mobility. At the same mobility in the yth direction, the permeation depth linearly decreases with anisotropic ratio as shown in Figure 10(a). The curve shape of permeation depth versus anisotropy ratio is similar for different mobility. In Figure 10(b), the permeation depth ratio decreases with the anisotropic ratio at the same mobility and increases with mobility at the same anisotropic ratio. The difference of permeation depth ratio is small regardless of mobility. The effect of permeation depth ratio can be fitted by Therefore, the correction coefficient in equation (16) can be expressed as Further, the permeation depth ratio is defined as where L α E =1 5,α ν =1 5 represents the permeation depth when α E = 1 5 and α ν = 1 5. Figure 11 presents the changes of permeation depth and permeation depth ratio with mechanical anisotropy at different mobility. Figure 11(a) shows that the permeation depth increases with mechanical anisotropy for either elastic modulus or Poisson's ratio at the same mobility. Larger mobility will have bigger permeation depth at the same mechanical anisotropy. In comparison with Poisson's ratio anisotropy, the elastic modulus anisotropy increases permeation depth more obviously. In Figure 11(b), the permeation depth ratio shows different tendencies at different mobility. As a whole, the elastic modulus anisotropy increases permeation depth ratio more significantly than the Poisson's ratio anisotropy. With the ratio ranging from 0.5 to 2, the permeation depth ratio changes from 0.99 to 1.01 for Poisson's ratio anisotropy, while it changes from 0.97 to 1.01 for elastic modulus anisotropy. The difference of permeation depth ratio is small regardless of mobility. The effects of permeation depth ratios in Figure 11 are fitted by Therefore, the correction coefficient in equation (16) can be obtained as f α E ,κ y = g α E f κ y = 1 004 − 0 112e −2 21α E −0 072 lg κ y − 0 365 , f α ν ,κ y = g α ν f κ y = 0 006α ν + 0 991 −0 072 lg κ y − 0 365 41 Equation (41) reflects the combined effect of mechanical anisotropy and mobility on the correction coefficient.

Conclusions
In this study, a numerical model was proposed to investigate the front motion of water permeation into an anisotropic porous medium. This model fully coupled the mechanical deformation, fluid flow, and moving water boundary in an anisotropic porous medium. This fully coupled model was verified by three sets of experimental data on concrete and shale samples. A revised formula with a correction coefficient was proposed for the quick estimation of the permeation depth of water front. Finally, parametric study was conducted to investigate the impacts of permeability, viscosity, permeability anisotropy, and mechanical anisotropy on the permeation depth of water front. The correction coefficient was analytically solved for radial flow pattern and porosity evolution and numerically evaluated for mobility, mechanical, and permeability 10 Geofluids anisotropy. Based on these preliminary studies, the following understandings and conclusions can be drawn.
(1) The proposed numerical model in this study has the capability to describe the front motion of water permeation into an anisotropic porous medium. The progress of permeation depth in both 1D and 2D water permeations can be well described by this model (2) Water mobility can well describe the combined effects of permeability and viscosity on the permeation depth of water front. The permeation depth increases with permeability, water mobility, and mechanical anisotropy but decreases with viscosity and permeability anisotropy. The impacts of mechanical anisotropy are smaller than those of water mobility and permeability anisotropy (3) A revised formula for the estimation of permeation depth was obtained through a correction coefficient in this study. The correction coefficient was analytically obtained for porosity evolution and radial flow pattern of water permeation. This correction coefficient was numerically calculated for the impacts of mobility, permeability anisotropy, and mechanical anisotropy, and curve fittings obtained can be used in engineering practice It is remarked that this model was established based on isothermal Darcy's flow, anisotropic and linearly elastic deformation. It did not consider the effects of cracking configurations, chemical reaction, non-Darcy flow, and temperature on water front motion. These are the topics in our future study. (c) Permeation depth ratio and fitting line is g α ν = 0 006α ν + 0 991 Figure 11: Changes of permeation depth and permeation depth ratio with mechanical anisotropic ratios.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.