Monte Carlo Assessment of the Impact of Oscillatory and Pulsating Boundary Conditions on the Flow Through Porous Media

Stress and water injection induce deformations and changes in pore pressure in the soil. The interaction between the mechanical deformations and the flow of water induces a change in porosity and permeability, which results in nonlinearity. To investigate this interaction and the impact of mechanical vibrations and pressure pulses on the flow rate through the pores of a porous medium under a pressure gradient, a poroelastic model is proposed. In this paper, a Galerkin finite element method is applied for solving the quasi-static Biot’s consolidation problem for poroelasticity, considering nonlinear permeability. Space discretisation using Taylor–Hood elements is considered, and the implicit Euler scheme for time stepping is used. Furthermore, Monte Carlo simulations are performed to quantify the impact of variation in the parameters on the model output. Numerical results show that pressure pulses and soil vibrations in the direction of the flow increase the amount of water that can be injected into a deformable fluid-saturated porous medium.


Introduction
Worldwide, management of freshwater resources has become critical. Climate scenarios predict extreme periods of drought and rainfall. Traditionally, during periods of heavy rainfall, the approach is to transport water quickly to surface waters and the sea in order to prevent flooding. However, a new approach is emerging: storing water in wet periods for use in dry periods. In particular, the storing of water underground has large benefits. A prerequisite for effective storing of rainwater in periods of extreme precipitation is that the water can be stored quickly. A new method to quickly infiltrate high volumes of fresh water has been discovered recently. We refer to this method as fast, high-volume infiltration (FHVI).
Another field in which FHVI constitutes a significant breakthrough is in construction, where this method was originally discovered. Building sites have to be pumped dry to enable construction. In the past, extracted water was often released into surface water. However, since the potential impact on the ecology is negative, new regulations prescribe that water should be returned to the ground. Conventionally, infiltration has to be applied at a certain distance from the pit, thereby affecting groundwater levels in a relatively large area. FHVI can be much closer to the pit, as the infiltrated water is rapidly transported away from the infiltration points. This means that the effect on groundwater levels is much smaller. However, according to preliminary research and findings, this infiltration method does not obey the classical Dupuit's law that is commonly used in hydrogeology, and it is currently impossible to predict its applicability.
In petroleum reservoirs, observations from the last 50 years suggest that seismic waves generated from earthquakes and passing trains may alter water and oil production. It has also been observed in some laboratory measurements and field applications that imposing harmonic signals into cores or reservoirs sometimes may induce higher fluid flow rates (Pan and Horne 2000). Further the fact that the pore pressure may undergo variations under the influence of seismic waves is well known to geotechnical engineers (Beresnev and Johnson 1994). In addition, laboratory experiments have shown that ultrasonic radiation can considerably increase the rate of flow of a liquid through a porous medium (Aarts and Ooms 1998). Furthermore, Davidson et al. (1999) performed experiments in a wide range of configurations, grain sizes, viscosities and flow factors, showing that high-amplitude pressure pulsing or mechanical excitation of a saturated porous medium under a pressure gradient increases the flow rate of the liquid along the direction of the flow gradient. Based on the experimental results, they concluded that flow rate enhancement occurred for all liquids, and for all the grain sizes that were tested.
Our aim is to investigate whether during FHVI, large injection rates induce an oscillatory or a pulsating force near the injection point and whether induced vibrations increase the amount of water that can be injected into the aquifer. In this publication, we tackle the second question by investigating the impact of soil vibrations and pressure pulses on the effective transmigration of water through the pores of the soil. For this purpose, a tube filled with a deformable fluid-saturated porous medium is simulated, into which water is injected. In the current paper, Biot's consolidation model for poroelasticity (Biot 1941(Biot , 1955) is used to determine the local deformation of the porous medium as a result of the injection of water. Biot's model is widely used in geomechanics, hydrogeology, petroleum engineering and biomechanics. Darcy's law (Hubbert 1957) is used in Biot's model to describe the fluid flow, while elasticity and poroelasticity of the porous medium determine the local deformations as a result of the vibrations. More precisely, we use in this paper a simplistic Hookean representation of the deformation of the soil.
The poroelasticity equations have been solved by Luo et al. (2015) using the finite volume method combined with a nonlinear multigrid method. In addition, stabilised finite difference methods using central differences on staggered grids are used by Gaspar et al. (2003Gaspar et al. ( , 2006 to solve Biot's model. However, the numerical solution of the two-dimensional poroelasticity equations is usually approached using finite element methods (Bause et al. 2017;Lewis and Schrefler 1978;Phillips and Wheeler 2007). In this study, a finite element method based on Taylor-Hood elements, with linear and quadratic basis functions, has been developed for solving the system of incompressible poroelasticity equations. This method is commonly used for flow problems modelled by (Navier-)Stokes equations. Furthermore, the fully coupled scheme was employed which involves solving the coupled governing equations of flow and geomechanics simultaneously at every time step. Another approach that is widely used in coupling the flow and the mechanics in porous media is the fixed stress split method (Both et al. 2017;Kim et al. 2009;Mikelić and Wheeler 2013). In this manuscript, we consider a nonlinear relation between the permeability and the dilatation. Subsequently, to quantify the impact of variation of model parameters such as Young's modulus, the oscillatory modes and the injection pressure pulses, we further present results from an uncertainty quantification. This uncertainty quantification is used to quantify the propagation of uncertainty in the input data. Such uncertainty quantifications have been applied in geomechanics (Luo 2017), where an uncertainty quantification is carried out by modelling the permeability as a stochastic field parameter.
The rest of this paper is organised as follows: Section 2 describes Biot's consolidation model. In Sect. 3, we formulate the numerical method. Here we derive the weak form of the partial differential equations and describe the Galerkin finite element approximations. Section 4 presents some of our numerical experiments and results. Lastly, in Sect. 5 we give our conclusions and make some suggestions for further work.

Model Equations
Sand and gravel layers (aquifers) are not rigid, but constitute an elastic matrix, if the deformations are very small. To be able to determine the local displacement of the skeleton of the porous medium, as well as the fluid flow through the pores, after injection of water, the model provided by Biot's theory of linear poroelasticity with single-phase flow is used (Biot 1941(Biot , 1955. In this model, flow in porous media is combined with mechanical deformations of the aquifer into which water is injected. Furthermore, Darcy's Law (Hubbert 1957) and infinitesimal strain theory (Bower 2010) are used to describe the fluid flow and the local displacements, respectively. Note that this is an approximation if the displacements and the strains are large.

Biot's Partial Differential Equations
The quasi-static Biot model for soil consolidation describes the time-dependent interaction between the displacement of the solid matrix and the pressure of the fluid. We assume the porous medium to be linearly elastic, homogeneous, isotropic and saturated by an incompressible Newtonian fluid. According to Biot's theory, the consolidation process satisfies the following system of equations (Aguilar et al. 2008;Wang 2000): (1) Darcy's law: Continuity equation: where σ is the effective stress tensor for the porous medium, p is the pore pressure, ρ is the density of water, g is the gravitational acceleration, λ and μ are the Lamé coefficients, u is the displacement vector of the porous medium, v is the percolation fluid velocity relative to the porous medium, κ is the permeability of the porous medium and η is the dynamic viscosity of the fluid.

Permeability and Porosity Relations
In the physical problem presented here, we will focus on the interaction between the mechanical deformations and the flow of water. Therefore, we consider the spatial dependency of the porosity and the permeability. We calculate the porosity using a procedure outlined by Tsai et al. (2006). Their derivation is based on the mass balance of solids in saturated porous media: where θ is the porosity and ρ s is the density of the solid skeleton. Assuming that ρ s is constant and that u is sufficiently smooth to interchange the order of differentiation with respect to time and space, we get ∂θ ∂t By Tsai et al. (2006), it is assumed that | ∂θ ∂t | |∇θ · ∂u ∂t |; herewith, they arrive at From this equation, we get with θ 0 the initial porosity. Note that the above relation differs from the one of Tsai et al. (2006), where in Eq. (7) they used a linearisation. By not applying this linearisation, we think that our approach is slightly more accurate. The permeability κ is determined using the Kozeny-Carman equation (Wang and Hsu 2009) where d s is the mean grain size of the soil. As a result of the dependency of the permeability on the mechanical deformations, problem (1)-(4) becomes nonlinear.

The Set-up of the Model
In this section, we will use Eqs.
(1)-(4) to describe the flow pattern of water in a tube filled with a poroelastic material, after the injection of water into the left end of the tube. The situation is as shown in Fig. 1a. We assume that the gravity-induced contribution to the flow of water is much smaller than the other contributions, which yields that the flow is axisymmetric, hence ∂ ∂θ (.) = 0 in whichθ is the azimuthal coordinate. Therefore, it is sufficient to determine the solution for a fixed azimuth (for example, the grey region in Fig. 1a). The computational domain Ω is thus a rectangular two-dimensional surface with cylindrical coordinates (x, r ), as depicted in Fig. 1b. Fig. 1 Sketch of the set-up for the tube problem: (left) physical problem and (right) numerical discretisation. Taking advantage of the symmetry of geometry and boundary conditions, only the grey region is discretised In order to solve this problem, Biot's consolidation model is applied on the computational domain Ω, with two spatial dimensions x = (x, r ) and with t denoting time: where˜ is the vector Laplacian. To complete the formulation of a well-posed problem, appropriate boundary and initial conditions are specified in Sects. 2.3.1 and 2.3.2.

Water Flow in a Vibrating Tube
To investigate the effect of vibrations on the water flow in the tube, we present two numerical experiments in this section. In these experiments, several ways of imposing vibrations are described, whereafter the effect on the volume flow rate at the right end of the tube is determined. In all problems, a tube of length L and initial radius R is considered. Furthermore, we assume that the casing of the tube is deformable, so that R = R(x, t) holds. The poroelastic material in the tube is assumed to be isotropic and homogeneous.
Effect of an Oscillating Casing of the Tube on the Flow In this problem, a tube is considered with a frictionless, impermeable casing on which transverse waves are imposed. Water is injected at a constant pressure into the soil through the left side surface (x = 0), while the right side surface (x = L) is kept at an ambient pressure at all times. Furthermore, filters are placed along the side surfaces to prevent that the grains exit the tube. More precisely, the boundary conditions for this problem are given as follows: where Γ = Γ 1 ∪ Γ 2 ∪ Γ 3 ∪ Γ 4 , t is the unit tangent vector at the boundary, n is the outer normal vector, u vib is a prescribed boundary displacement due to the vibrations and p pump is a prescribed boundary pore pressure due to the injection of water. Figure 1b shows the definition of the boundary segments. Note that the boundary conditions on boundary segment Γ 3 are required by the definition of symmetry. The variational inequality in condition (11g) accounts for the fact that the grains cannot exit the tube through boundary segment Γ 4 . More specifically, condition (11g) states: u·n ≤ 0 and (σ n)·n = 0 or u·n = 0. This boundary condition could also be used on boundary segment Γ 2 . However, in this case it is possible that (σ n) · n = 0 on both boundary segments Γ 2 and Γ 4 . Then, there is no Dirichlet boundary condition for the horizontal displacement (in the x-direction). This leads to a degenerate elliptic operator for the displacement u, which could make the problem ill-posed. Initially, the following condition is fulfilled: As mentioned earlier, transverse waves are imposed for the boundary displacement u vib , represented by with γ the amplitude of the wave, λ w the wavelength and v the phase velocity. Note that for v < 0 the wave is travelling to the left, while for v > 0 the wave is travelling to the right.
Effect of a Vibrating Imposed Load on the Flow While in the previous section the vibrations were imposed as an oscillating casing of the tube, in this section the effect of an oscillating load applied on the casing is analysed. Accordingly, the boundary condition for the mechanical deformation on Γ 1 becomes with σ vib is an oscillating vertical load. Similar to the previous section, transverse waves are used for the oscillating load On t = 0, the initial condition (12) is fulfilled.

Pulsed Injection
Instead of applying mechanical vibrations on the displacement, we can also investigate the effect of a pulsed injection of water into the left end of the tube. In this case, the prescribed boundary pore pressure p pump caused by the injection of water will have a pulsating behaviour rather than being constant. Hence, for the boundary conditions for the mechanical deformation on Γ 1 and for the water pressure on Γ 4 holds with p pump (t) represented by the Heaviside step function A rectangular pulse wave with period T p and pulse time τ can now be defined as where p max is the maximum injection pressure and N p is the number of periods. Furthermore, the initial condition (12) is fulfilled.

Weak Formulation
To present the variational formulation of these problems, we first introduce the appropriate function spaces. Let L 2 (Ω) be the Hilbert space of square integrable scalar-valued functions f We further introduce the function space Q = {q ∈ H 1 (Ω) : q = 0 on Γ 2 and q = p pump on Γ 4 } for all the problems that we consider. Subsequently, we use the function spaces W = {w ∈ (H 1 (Ω)) 2 : w·n = u vib on Γ 1 and w· n = 0 on Γ 2 } and W 0 = {w ∈ (H 1 (Ω)) 2 : w · n = 0 on Γ 1 ∪ Γ 2 } for the problems with boundary conditions as stated in Eqs. (11) and (16), respectively. For the problem with boundary conditions (14), Furthermore, we consider the bilinear forms (Prokharau and Vermolen 2009) The variational formulation in cylinder coordinates (x, r ) for problem (10) with boundary and initial conditions (11)-(12) and also for problem (10) with initial and boundary conditions (12) and (16) consists of the following, using the notationu = ∂u ∂t : For each t > 0, find with the initial condition u(0) = 0, and where The variational formulation for problem (10) in cylinder coordinates (x, r ) with initial and boundary conditions (12) and (14) consists of the following: with the initial condition u(0) = 0, and where

Finite Element Discretisation
Problems (21)- (24) and (25)-(27) are solved by applying the finite element method, with triangular Taylor-Hood elements (Braess 2001;Van Kan et al. 2008;Segal 2012). Let P k h ⊂ H 1 (Ω) be a function space of piecewise polynomials on Ω of degree k. Hence, we define finite element approximations for W and Q as (Aguilar et al. 2008;Prokharau and Vermolen 2009). Subsequently, we approximate the functions u(t) and in which the Dirichlet boundary conditions are imposed. Then, the semi-discrete Galerkin approximation of problem (21)-(24) is defined as follows: and for t = 0: u h (0) = 0. Simultaneously, discretisation in time is applied using the backward Euler method. Let t be the time step size and define a time grid {t m = m t : m ∈ N}, then the discrete Galerkin scheme of (29)-(30) is formulated as follows: while for m = 0: u 0 h = 0. The discrete Galerkin scheme for problem (25)-(27) is derived similarly. These discrete Galerkin schemes are solved using triangular Taylor-Hood elements. The displacements are spatially approximated by quadratic basis functions, whereas continuous piecewise linear approximation is used for the pressure field. We remark that the Taylor-Hood elements are suitable as a stable approach for this problem. However, spurious oscillations are diminished but not completely removed for small time steps. To fully remove the non-physical oscillations, one may use the stabilisation techniques as considered by Aguilar et al. (2008). The numerical investigations are carried out using the matrix-based software package MATLAB (version R2011b). At each time step, we solve Eqs. (31)-(32) as a fully coupled system, where we use the permeability from the previous time step. After having obtained the numerical approximations for u and p, we update the porosity using Eq. (8). Subsequently, the Kozeny-Carman relation (9) is used to calculate the permeability. The new value for the permeability is then used for the next time step. An iterative method is not used in this approach because of efficiency and since no instability was observed in our results. In order to use Eq. (8), we determine the dilatation using the numerical solution u h = (u x,h , u r,h ) at time t m . The spatial derivatives in the dilatation are then computed by applying the finite element method. Firstly, we introduce the functions d ∈ L 2 (Ω) and we define the finite element approximation as Hence, the discrete Galerkin scheme is given by The discrete Galerkin scheme for the function ω m 2 is derived similarly. In these finite element schemes, the spatial derivatives are approximated by quadratic basis functions. For the integrals in the element matrices and the element vectors, exact integration is used. Regarding accuracy, our numerical experiments showed that this strategy produced sufficiently reliable results. We note that improvements of the accuracy can be obtained using gradient recovery techniques which yield superconverging behaviour (Zienkiewicz and Zhu 1992). The aim of this study is to investigate the impact of oscillatory and pulsating boundary conditions on the volume flow rate at the right end of the tube. We compute the volume flow rate at the right boundary segment, using the velocity field as described by Darcy's law (3). To compute the velocity field for post-processing issues, the finite element method is applied analogously to the computation of the derivatives of the displacements, combined with piecewise linear approximation and exact integration.

Numerical Simulations
In this section, we discuss the solution results for the discrete Galerkin approximation of the quasi-two-dimensional problems that are presented in Sect. 2.3. The simulation domain is a rectangle with length 1.0 m and width 10 cm (see Fig. 1b). The domain is discretised using a 101 by 11 regular triangular grid, which provides sufficient resolution according to a mesh refinement study. The chosen values for the material properties of the porous medium are given in Table 1, where λ and μ are related to the Young's modulus E and the Poisson's ratio ν by (Aguilar et al. 2008 1+ν) . The values in Table 1 are chosen based on discussions with experts from engineering and consultancy company Fugro GeoServices B.V. and on the literature (Van Wijngaarden 2015).

The Impact of an Oscillating Casing of the Tube on the Water Flow
In order to obtain some insight into the impact of an oscillating casing of the tube on the water flow, we present an overview of the simulation results in Figs. 2 and 3. In this simulation, water is injected into the soil at a constant pump pressure equal to 0.5 bar. As a consistency   Table 1. a Numerical solution for the pressure. b Numerical solution for the fluid velocity. c Numerical solution for the permeability. d Positions of the mesh points after subjecting them to the calculated displacement vector u the right, as a result of the force exerted on the grains by the injected water. As a result of this small shift of the grains and the assumption that the grains cannot exit the tube, we expect a higher grain density near the right end. Consequently, the permeability will linearly decrease towards the right end of the tube, as depicted in Fig. 2c. In Fig. 3, the numerical solutions are shown for a test case of problem (31)-(32) with vibrations. In this test case, a transverse wave (13) travelling to the right is chosen as prescribed boundary displacement u vib , with γ = 1 cm, λ w = 1 m and v = 1 m/s. In contrast to the pressure shown in Fig. 2a, the numerical solution for the pressure in the problem with vibrations is no longer linear, but shows an oscillatory behaviour, as depicted in Fig. 3a. The vibrations also provide an oscillatory profile in the permeability, as shown in Fig. 3c. In this figure, we can see that the permeability decreases when the grains are pressed together by the vibration, while it increases when the grains are pulled apart. The simulation results in Figs. 2 and 3 show an impact of the vibrations, imposed on the casing, on the water flow. However, by only looking at these results, the impact of vibrations on the amount of water that flows through the tube stays unmeasurable. For this reason, the impact of the vibrations and pulses on the water flow is defined in this publication as the  Table 1 impact on the volume flow rate Q at the right end of the tube. In Fig. 4, two graphs are presented that sketch the behaviour of the volume flow rate over time at the right end of the tube. In these simulations, the aforementioned test cases are used. For the test case without vibrations, the time average of the volume flow rate Q over the time interval (0, 5] is equal to 6.86 × 10 −5 m 3 /s. For the test case with imposed vibrations, the time average of the volume flow rate Q over the time interval (0, 5] is equal to 9.69 × 10 −4 m 3 /s. Thus, the percentage change Q % of the time average of the volume flow rate as result of the imposed vibration is 1311.7%. Based on these test cases, we can conclude that the volume flow rate at the right end of the tube increases as a result of the imposed vibrations on the casing, as depicted in Fig. 4. We finally note that the area enclosed by the Q-curve and the t-axis represents the total amount of water that flows out of the domain over a certain period.

Monoparametric Variation
In this section, we will investigate the impact of the transverse waves (13) Fig. 5, the time step size t is determined using the formula t = 1 8 f , in which f is the frequency computed by: f = |v| λ w . In case variation is applied to the values of λ w or v, the maximum value of f is used to determine the time step size. Figure 5a indicates that an increase in the amplitude of the imposed wave, with fixed wavelength and phase velocity, results in an increase in the volume flow rate. However, an increase in the wavelength, with fixed amplitude and phase velocity, leads to a decrease in the volume flow rate, as shown in Fig. 5b. Figure 5c  in the phase velocity magnitude |v|, with fixed amplitude and wavelength, leads to a larger impact on the volume flow rate. Furthermore, assigning positive values to the phase velocity v (corresponding with transverse waves travelling to the right) results in positive volume flow rate profiles. Assigning negative values to the phase velocity v (corresponding with transverse waves travelling to the left) produces negative volume flow rates. As expected, the water flow, which is directed to the right, is stimulated by waves travelling in the same direction. However, waves travelling in the opposite direction counteract the flow, resulting in negative volume flow rates. In fact, the negative volume flow rates are a result of the force applied by the oppositely directed waves that is larger than the pump pressure. At a higher pump pressure, the effect of these waves on the volume flow rate is smaller, as illustrated in Fig. 6.

Monte Carlo Method
In the previous section, the contributions of the variations in the values of the vibration characteristics to the volume flow rate were quantified by applying a monoparametric variation where the values of the parameters are varied one by one, within ranges of possible values. As we choose to fix a number of parameter values each time, we are not able to draw any conclusions from this monoparametric variation. In this section, Monte Carlo method is applied to the values of the vibration characteristics using samples of uniform distributions with chosen boundaries. In Fig. 7, the time average of the volume flow rate is depicted after applying Monte Carlo simulations to the values of the vibration characteristics γ , λ w and v. For the simulations, 300 samples from the following uniform distributions are generated:  Table 1 γ ∼ U (0.001, 0.01), λ w ∼ U (1/4, 1), v ∼ U (−2, 2). Figure 7a shows that an increase in the amplitude of the imposed wave leads to a larger impact on the time average of the volume flow rate Q. However, an increase in the wavelength results in a smaller impact on Q, as shown in Fig. 7b. In both figures, we observe a mirroring across a line near the horizontal axis. Values of Q above the mirror line correspond with positive values of the phase velocity, while values of Q bellow the mirror line correspond with negative values of v, as can be concluded from Fig. 7c. Furthermore, Fig. 7c shows that an increasing phase velocity magnitude |v| leads to a larger impact on Q. In addition, we observe in Fig. 7c that for v < 0 some of the values of Q are positive, and these values correspond with small values of the amplitude. In this case, the pump pressure is larger than the force applied by the oppositely directed waves, which results in positive volume flow rates. Note that these results are consistent with the results from Sect. 4.1.1. From the formula for the frequency and Fig. 7, we can conclude that large amplitudes and high frequencies lead to high-volume flow rates for v > 0. In Table 2, the Pearson correlation coefficients are given together with the associated p-values. From this table, we can conclude that the vibration characteristics γ and v have the most impact on the time average of the volume flow rate.
To measure the real impact of the transverse waves (13) on the water flow, the percentage change Q % of the time average of the volume flow rates as result of the imposed vibrations is determined by the formula Q % = Q−Q 0 Q 0 × 100, where Q 0 is the time average of the volume flow rate in the test case without vibrations. The percentage change Q % of the time average Q that is computed after applying the above-described Monte Carlo method (see Fig. 7), is depicted as function of γ in Fig. 8. Figure 8 shows that, for v > 0, the volume flow rate at the right end of the tube increases as a result of the imposed vibrations on the casing. For vibrations with a large amplitude and a high frequency, the time average of the volume flow rate can become as large as 42 times the time average of the volume flow rate in the test case without vibrations. The smallest percentage change in the volume flow rate as cause of the vibrations is equal to 6.0%. On the other hand, for v < 0, all vibrations lead to a negative percentage change Q % , even the vibrations with small values of the amplitude γ . Given the probability space (Ω p , F , P), with sample space Ω p (which is the set of all possible outcomes), the set of events F (where each   event is a set with zero or more outcomes), and the probability P, with P : F → [0, 1], we can compute the cumulative distribution function. In Fig. 9, the histogram and the cumulative distribution function of the percentage change Q % are presented. Since 51% of the sampled values of v are negative, we see in Fig. 9 that P(Q % ≤ 0) ≈ 0.51. Furthermore, using the probability space we can compute the probabilities that the percentage change is greater than 100%: P(Q % ≥ 100) ≈ 0.43 and P(Q % ≥ 100|v > 0) ≈ 0.89.

The Impact of a Vibrating Load on the Water Flow
In this section, the impact of a vibrating load, which is applied on the casing of the tube, on the water flow is investigated. This numerical experiment makes it possible to analyse the contributions of the variations in the values of the porous medium properties to the volume flow rate. These contributions are quantified by assigning a range of possible values to the parameters: E, ν, θ 0 and d s . For this purpose, four types of soil are distinguished: clay, silt, sand and gravel. In the literature, there is a large consensus that the Kozeny-Carman equation (9) applies to sands but not to clays (Chapuis and Aubertin 2003). Therefore, this experiment is only applied to sand and gravel.

Monte Carlo Method for Sand and Gravel
Monte Carlo method is applied to the values of the material properties of sand and gravel, using samples of uniform distributions with boundaries found in the literature ("Geotechnical parameters" 2017; "Soil properties" 2017). As the values within these boundaries are all equally likely to occur, we have chosen to use uniform distributions instead of another frequently used distributions like the log-normal distribution. In Figs Table 1 while for gravel, 300 samples from these uniform distributions are generated: E ∼ U (80 × 10 6 , 160 × 10 6 ), ν ∼ U (0.30, 0.40), θ 0 ∼ U (0.23, 0.38), d s ∼ U (2 × 10 −3 , 50 × 10 −3 ).  Figures 10 and 11 show that water flows faster through porous media with large mean grain sizes or high initial porosities, after imposing a vibrating load on the casing. On the other hand, these figures show that the volume flow rate is invariant under variation in the values of Young's modulus and Poisson's ratio. This can also be concluded from the values of the Pearson correlation coefficients given in Table 3.

The Impact of a Pulsed Injection on the Water Flow
Based on over 160 laboratory tests, Dusseault (1999) demonstrated the beneficial effects of pressure pulsing on the water flow in porous media. To theoretically examine his findings, we investigated the impact of a pulsed injection of water, into the left end of the tube, on the water flow. The results of this research are presented in this section. In addition, the contributions of the variations in the values of the pulse wave characteristics to the volume flow rate are analysed. These contributions are examined by applying Monte Carlo method to the pulse wave period T p and the relative pulse timeτ , defined asτ = τ T p , with τ the pulse time (see Formula (18)). Subsequently, in order to be able to draw reliable conclusions about the impact of pressure pulsing on the volume flow rate, we compare the volume flow rate caused by the pulsed injection with the volume flow rate by a constant pump pressure p pump . In this comparison, the percentage change Q % of the time average of the volume flow rate Q is used, determined by Q % = Q−Q C Q C × 100, where Q C is the time average of the volume flow rate caused by a constant pump pressure. In each simulation, the constant pump pressure is chosen equal to the total pump pressure over time by pulsed injection. Hence, for a particular relative pulse timeτ , the constant pump pressure is computed by p pump =τ p max . In Fig. 12, the percentage change in the time average of the volume flow rate at the right end of the tube is depicted after applying Monte Carlo simulations to the values of the pulse wave characteristics T p andτ . For the simulation, 300 samples from the following uniform distributions are generated: In the generations of the simulation results presented in Fig. 12, the time step size t is determined by t = T p 20 . As variation is applied to the values of T p , the minimum value of t is used as time step size. Figure 12 shows that pressure pulsing with small relative pulse timesτ leads to a major increase in the volume flow rate, while it increases slightly by increasing the pulse period T p . This can also be confirmed by the Pearson correlation coefficients and the p-values given in Table 4. In Fig. 13, the histogram and the cumulative distribution function of the percentage change Q % are depicted. Since P(Q % ≤ 0) = 0, we can conclude from Figs. 12 and 13 that pulsed injection has a beneficial effect on the water flow in porous media.  Table 1. a Q % as function of T p . b Q % as function ofτ

Conclusions and Discussion
In this study, the poroelasticity system with nonlinear permeability is solved using the Galerkin finite element method based on Taylor-Hood elements, combined with a backward Euler time integration. The study contains simulations with oscillatory force boundary conditions as well as pressure pulses. Furthermore, to quantify the impact of variation of model parameters such as Young's modulus, the oscillatory modes and the injection pressure pulses, a probabilistic approach is carried out.
To begin with, soil vibrations are applied on the casing of a tube as oscillatory displacement boundary condition. Numerical results showed that large amplitudes and high frequencies of the imposed mechanical vibrations lead to high-volume flow rates for positive values of the phase velocity, corresponding with transverse waves travelling to the right. Therefore, the volume flow rate at the right end of the tube increased as a result of the imposed mechanical vibrations. On the other hand, for negative values of the phase velocity, corresponding with waves travelling to the left, the vibrations lead to a decrease in the volume flow rate at the right end of the tube. As expected, the water flow, which is directed to the right, is stimulated by waves travelling in the same direction, while waves travelling in the opposite direction counteract the flow, resulting in negative volume flow rates in case the force applied by the oppositely directed waves is larger than the pump pressure.
Subsequently, applying an oscillating load on the casing of the tube showed that water flows faster through porous media with large grain sizes and/or high initial porosities. On the other hand, variation in the values of Young's modulus and Poisson's ratio indicated that these parameters do not have a large impact on the volume flow rate.
Numerical simulations of pressure pulsing pointed out that injection pulses with small relative pulse times lead to a major increase in the volume flow rate, while an increasing pulse period results in a slight increase in the flow rate. Most importantly, we can conclude that pulsed injection has a beneficial effect on the water flow in porous media.
In the current paper, we use a "brute-force" Monte Carlo simulation procedure. Doing simulations with thousands of samples can reduce the Monte Carlo error. However, as in our case each sample simulation takes about two hours, we instead adopted 300 samples. In recent studies (Cliffe et al. 2011), multilevel Monte Carlo methods (MLMC) have been applied to various systems, such as groundwater flow. These MLMC methods are characterised by the advantage that relatively few simulations are needed at high mesh resolutions, whereas one performs large numbers of simulations at lower resolutions. The MLMC methods are therefore thought to be a suitable candidate for future applications.
In conclusion, pressure pulses and soil vibrations in the direction of the flow increase the amount of water that can be injected into a tube filled with a deformable fluid-saturated porous medium. However, to elucidate the underlying mechanisms of FHVI, the injection of water into the aquifer should be simulated. In addition, further research should be carried out to examine the validity and applicability of the Kozeny-Carman equation especially for the problems in which vibrations can lead to very large gradients of the porosity and permeability in the soil. Furthermore, it is important to note that the model we are currently using is based on the assumption that the displacements and the strains are small. For this reason, we used Hooke's law and the assumption that the strain tensor is only determined by the deformation gradient (and its transpose), while higher-order terms are neglected. However, as the mechanical vibrations can lead to large strains and displacements, it is probably more realistic to use a morphoelastic model or another plasticity model for soil. The treatment of large deformations would imply nonlinear terms. For the morphoelastic model, for example, an additional nonlinear time-dependent equation would have to be solved, which would increase the computational complexity. Since the purpose of this study is to investigate how vibrations can influence the fluid flow using uncertainty quantification, these nonlinear terms are neglected and the model is kept simple and cheap from a computational point of view.