Multiphysics modelling of the Mandel-Cryer effect

In porous medium studies the Mandel-Cryer effect is known, describing non-monotonic pore-water pressure evolution in response to loading or to changed stress conditions. In a 2D poro-elastic model we couple the pore water hydraulics with mechanics (HM). The Mandel-Cryer effect is identified in parts of the model region that are far from the drainage boundary. At parts of the loaded boundary an even more complex pressure evolution is revealed. Variations of the Biot-parameter as the coupling parameter clearly indicate the relevance of the two-way coupling between the involved physical regimes. Hence the Mandel-Cryer effect is a typical result of multi-


INTRODUCTION
In problems of poro-elasticity the so-called Mandel-Cryer effect describes the nonmonotonic response of pressure due to external loading.A sudden change of pressure on an open system usually is expected to result in a monotonic reaction: the initial pressure will gradually decay towards a new equilibrium state.However, in a poro-elastic system an increase of pressure can be observed in parts of the system for a certain time, before the overall decrease becomes dominant.
Mandel [1] was the first, who described such behaviour.Mandel's problem consists of an infinitely long rectangular specimen sandwiched at the top and the bottom by two rigid plates.The lateral sides are free from normal and shear stress, and pore pressure.At t = 0 a force is applied to the rigid plates and a uniform pore pressure appears, according to the Skempton effect [2].With progressing time the pressure near the side edges will dissipate because of drainage.But initially the solution, derived by Mandel, shows an increase of pressure.The transferring of compressive total stress works as a pore pressure generation mechanism such that the pressure in the centre region continues to rise after its initial creation.In the centre between the drained ends the pressure increase is most pronounced and also prevails the longest time.
Altogether the pore pressure response hence is non-monotonic, a characteristic not observable in a simple diffusion phenomenon such as that modelled by the Terzaghi theory [3].It has to be noticed that the Mandel set-up is 1D, as details pressure and deformation development are described only in the single space dimension, which extends from the centre of the specimen to the drained ends.In our numerical model we will elaborate on a 2D set-up, in which the processes in the transverse (vertical) direction are included also.After Mandel, Cryer [4] observed a similar behaviour at the centre of a sphere consolidating under hydrostatic pressure.The sphere is open at the surface, so that the initial sudden pressure increase will be compensated by fluid drainage.Cryer developed an analytical solution for this situation, which in the interior of the sphere shows the same non-monotonic development of pressure, as described by Mandel before.
This type of non-monotonic pressure response is usually referred to as the Mandel-Cryer effect.It is a distinctive feature of the multi-physical coupling, in contrast to the traditional uncoupled Terzaghi theory [3].The physical phenomenon has been confirmed in the laboratory [5,6].Abousleiman et al. [7] extended the Mandel problem to include transverse isotropy as well as compressibility of pore fluid and porous medium (see also: Cui et al. [8]).
Finite element modelling is a convenient numerical technique for the simulation of the Mandel-Cryer effect.They have been applied for the Mandel problem [8,11,12,17], as well as for the Cryer problem [18,19].Including thermal effects additionally Selvadurai & Suvorov [20] extend the Cryer setting to a coupled THM problem, also using the code COMSOL Multiphysics [21].Using finite elements Jha [22] observed and described the Mandel-Cryer effect in a constellation of groundwater pumping and injection wells.Here, we focus on a modified version of the Mandel set-up in 2D that was proposed by Jha & Juanes [15].

HYDRAULIC-MECHANICAL COUPLING
The Mandel-Cryer effect results from hydraulic-mechanical coupling (HM).The hydraulic regime in question is porous medium flow according to Darcy's Law, determined by hydraulic pressure in the pore space.The pressure however is also affected by deformations of the elastic material, described by Biot's theory.Moreover the influence is vice versa, as the geomechanics depends on pressure, too.Thus there is a two-way coupling.In the sequel the differential equations for the single processes are described first and the coupling terms are introduced thereafter.
The stress regime the porous medium and the corresponding deformations are described by the mechanical equations Eqs ( 1) and ( 2) is a system of differential equations in which the elements of the deformation vector u are the dependent variables.Further variables are: stress tensor σ, volume force vector   and stiffness tensor C. The strain tensor ε is defined by Equation ( 2) describes the stress-strain relationship in general.In the application below we deal with elastic media only.Subscript 0 denotes initial or reference values.
The hydraulic conditions in porous media are described by: Equation ( 4) is a differential equation for pore pressure p as dependent variable.Other variables and parameters are: Darcy velocity vector q, fluid density ρ, permeability k, fluid dynamic viscosity μ, acceleration due to gravity g, storage parameter S, and fluid source/sinkterm Q. Equation ( 4) is a formulation of the mass balance and equation ( 5) is Darcy's Law for flow in porous media.The equations ( 4) and ( 5) describe the diffusion of pressure, where the diffusivity D is given by: The pore pressure also affects the mechanical state of the system.In the mathematical description this is taken into account by an extension in the stress-strain equations (2): with Biot constant α [23,24].For high Biot constants (≅ 1) there is a strong coupling, for low α (≅ 0) there is a weak coupling.In order to consider anisotropy direction-dependent Biot constants can be introduced in the last term on the right hand side of the equation [7].Cui et al. [8] use a tensor of Biot-constants as a further generalization.The backward link from mechanics to hydraulics is given by an additional term in the fluid equation ( 4): with volumetric strain ε v .The term also depends on the Biot constant α.

MANDEL-CRYER MODEL SET-UP
A model was set up that illustrates the Mandel-Cryer effect in the coupled poro-elastic approach described above.The original sandwiched specimen used by Mandel [1], showing the non-monotonic behaviour of pore pressure, was modified.While Mandel's setting becomes basically a 1D problem if transversal effects are neglected, the latter are considered in our 2D approach.In the conceptual model we follow Jha & Juanes [15].While in the original Mandel set-up only vertical loading was assumed, the Mandel-Cryer model, used here, is driven by lateral loading from the side.Using overburden and sideburden Kim et al. [16] propose another very similar 2D modification of the original Mandel set-up.
In the model a specimen is exhibited to a uniform constant compressive pressure of σ 0 =1 MPa from the left.The right and bottom boundaries are fixed in the normal directions, but allow movements in tangential directions ('roller' conditions).Concerning flow, the top is a drained boundary with constant pore pressure p=0, while the other three boundaries are no-flow boundaries.The situation is sketched in Figure 1 (displacement vector u = (u,v)).For modelling we use the software COMSOL Multiphysics [21], a versatile and flexible code for Finite Element solutions of coupled partial differential equations.The software can be applied to all kinds of (multi)-physical settings.It is equipped with a graphical user interface that allows easy handling and coupling of different physics modes.There are several toolboxes that extend the core program.For this study we utilized the 'porous medium and subsurface' toolbox and the 'poroelastic' mode.
With the 2D deformation vector u and pore pressure p there are three dependent variables in the coupled system of equations.In the finite element formulation all dependent variables are discretized by quadratic element functions on a triangular mesh.The finite element mesh is refined at the top boundary.Altogether we used 72139 quadratic elements, 5057 edge elements, with an average element quality 0.95.

RESULTS
The model was simulated for the time interval from 0 to 1400 s. Figure 2 shows the deformation of the model, i.e. its shrinking in horizontal direction and its extension in vertical direction after 1000 s for the maximum coupling case, which we use here as a reference.The length unit on the colorbar is [m].For illustration purposes the displacement is exaggerated.
Figure 3 depicts the development of pressure with time for the reference case with maximal coupling, i.e. α=1.The pore pressure at the right and lower boundaries is plotted for different  The Mandel-Cryer effect, the early times pressure increase, can clearly be recognized at parts of both boundaries.At the lower part of the vertical boundary and most parts of the bottom boundary, the pressure increases before beginning to decrease.At t=0, an undrained pressure is generated.The specimen is deformed: due to the boundary conditions it shrinks in horizontal direction and extends in vertical direction upwards.With progressing time the pressure near the top boundary decreases because of fluid drainage.As the hydraulic diffusivity is small, the effect of drainage is not observed immediately near all no-flux boundaries.This results into load transfer of compressive total stress towards the bottom boundary, in response to which the pressure there continues to rise above its undrained value.
Near the left boundary, where the specimen is loaded, the situation is more even more complex.At early times we observe a pressure decrease due to drainage in the direct vicinity of the loading.The system then adjusts to the increased pressure regime at the bottom, that was described, which leads to a pressure increase in a second phase, before the decrease in the final phase.The details of pressure evolution at three positions at the bottom boundary, shown in Figure 4, illustrate the complex behaviour.In the long term all excess pressure vanishes and a uniform stress regime returns.Hence, the pressure evolution at points away from the drained boundary is non-monotonic, a phenomenon not observed in a purely diffusive process such as that modelled by the Terzaghi theory, where the pressure is uncoupled from the solid deformation.
In order to study the coupling effect we ran several model runs with different values for the Biot parameter.The coupling parameter was changed from α=0 (no coupling) to α=1 (maximal coupling).Figures 5 shows the evolution of pore pressure in the lower right corner of the model region for all α.
For the uncoupled case the pressure remains constant throughout: the pressure increase at one side of the model is swiftly compensated by drainage at the upper boundary.For coupled regimes the non-monotonic increase can be observed at the lower right corner in all cases.p is increasing before it is declining.For α=0.2 the pressure increase is hardly visible.However, the non-monotonic behaviour can be clearly identified as a coupling effect.The maximum pressure is reached at t=100 s for α =0.2, at t=178 s for α =0.4, at t=245 s for α =0.6, at t=287 s for α =0.8 and at t=355 s for α =1.The coupling strength obviously extends the time of pressure build-up.While the peak size increases less for strong coupling, the relation between peak time and Biot-parameter is nearly linear.
The representation at logarithmic time scale in Fig. 5 reveals that after loading there is a time period in which the pressure remains quasi constant, before a rise becomes apparent.Thus for the lower right corner of the model in fact three different phases have to be distinguished: first an initial period with constant pressure, a second with rising pressures and a third with decreasing pressures.
Fig. 6 displays the maximum deformation at the right boundary for different coupling parameters.We choose the same time, here: 1400 s.As expected the deformation increases with the coupling.However the relation between deformation and α is nonlinear.For weakly coupled cases the deformation changes slightly with α, while for strongly coupled cases the change of deformation is much bigger.

CONCLUSIONS
The Mandel-Cryer effect, the initial pressure increase as response to an external load, is studied as a multi-physics problem concerning a 2-way coupled regime of hydraulics and mechanics.We studied the effect in an extended 2D set-up of the original Mandel constellation, which allows the study of pressure development in transverse direction in addition.
The non-monotonic behaviour of pressure can be observed in the lower half of the system only.The extension of the region with Mandel-Cryer effect surely depends on the parameter values.Moreover near to the loading boundary we observe an even more complex behaviour: in the lower part the pressure reacts in three phases: decrease, increase and final decrease.The development can be explained by the different physical regimes, which become dominant at different time periods.

Figure 1 :
Figure 1: Sketch of the 2D modified Mandel set-up instants (the legend shows time in [s]).The y-coordinate increases from bottom to the top, the x-axis from left to the right.Pressure is non-dimensionalised with the applied compressive stress.

Figure 2 :Figure 3 :
Figure 2: Total displacement for maximal coupling after 1000 s; for the visualization the deformations at left and upper boundaries are exaggerated

Figure 4 :
Figure 4: Evolution of (non-dimensional) pressure at selected positions at the bottom boundary; for x=0 the pressure curve with a minimum and a maximum shows more complex development

Figure 5 :
Figure 5: Pore pressure evolution at the lower-right model corner for coupling parameters α=0:0.2:1;upper sub-plot for linear time scale, lower sub-plot for logarithmic time scale

Table 1 :
Parameters for the 2D modified Mandel model