Skip to main content


Front. Phys., 18 March 2022
Sec. Interdisciplinary Physics
Volume 10 - 2022 |

Drainage Instabilities in Granular Materials: A New Biaxial Apparatus for Fluid Fingering and Solid Remodeling Detection

www.frontiersin.orgRana Al Nemer* www.frontiersin.orgGiulio Sciarra www.frontiersin.orgJulien Réthoré
  • Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR, Nantes, France

In this study, the effect of fluid fingering on the solid remodeling of a granular material during the drainage phenomenon is investigated. A new biaxial apparatus, endowed with two transparent windows and adapted to unsaturated soils, is used to capture the effects of hydraulic instabilities on the mechanical response, by means of high-resolution cameras. A specimen of (40 × 50 × 11) mm3, of Fontainebleau sand NE 34 initially saturated by water, is connected to a pressure-controlled gas source to inject the gas into the sample. During the injection phase, fluid instabilities are detected and filmed. Using imaging techniques, the grain remodeling and strain localization due to the two-phase fluid flow are measured.

1 Introduction

Fluid fingering and strain localization are two kinds of instabilities possibly occurring in granular media when a fluid front percolates the porous network or when a perturbation of an initial state of equilibrium is induced by a non-vanishing deviatoric state of stress, respectively. Such kind of phenomenon is ubiquitous in geomaterials. Fluid fingering instabilities are involved, for instance, in the long-term fate of stored underground hydrocarbons or sequestrated CO2 (see [1]), leakage of non-aqueous phase liquids in unsaturated granular media (see [2]) or drying of initially saturated soil (see [3]), strain localization instabilities in fault reactivation (see, e.g., [4, 5]), and shear band or fracture nucleation (see, e.g., [68]).

Since the eighties, two separate scientific communities have worked on these topics, the first being mostly interested in understanding the flow of multiple immiscible phases through a rigid porous skeleton, and the second in characterizing the strain and stress localization in granular materials in dry or fully saturated conditions, under biaxial or triaxial loading.

In what concerns multi-phase flow through porous media, the forementioned fluid instabilities have been investigated considering a network of capillary ducts [9] and more recently [10], Hele-Shaw cells [11], or real porous rocks [12,13], attempting at finding out a relation between pattern formation at the scale of centimeter-scale samples and phenomena occurring at the scale of the pores, such as piston-type motions, snap-offs, and Haines jumps. Indeed, different pore-scale events show whether the multi-phase flow occurs in a drainage or an imbibition process. The focus of this study is on drainage processes, when a wetting phase (WP) initially saturating a granular medium is displaced by the injection of a non-wetting phase (NWP). At the pore scale, the invasion of the NWP usually corresponds to an alternate sequence of piston-type displacements and Haines jump events when the NWP moves from pores to throats and from throats to pores, respectively. The characteristic time scales of these two processes are very different; as a matter of fact, pore filling due to Haines jumps abruptly occurs, as a kind of pore-scale instability, with a sharp drop of the local capillary pressure caused by the change in the curvature of the non-wetting/wetting interface. The fast invasion of the pore space by the NWP generally induces a non-wetting fluid flow toward the pore from adjacent throats, hence sucking wetting fluid into these last ones. These pore-scale instability events are clearly driven by capillary forces and therefore happen only when capillary forces exceed the viscous force. This means, referring to the Lenormand phase diagram [9], that Haines jumps mostly occur within the capillary fingering domain [14, 15]. On the other hand, the regime of viscous fingering, with viscous forces higher than capillary forces, typically implies, at the microscopic scale, a piston-type invasion, driven by high values of the NWP flow and, consequently, of the pressure gradient, this last being responsible for the absence of the stick–slip behavior previously described.

For what concerns strain localization and crack/fault (re-)opening in geomaterials, the focus has been on the quantitative characterization of such instabilities (extent, thickness, and void index evolution, see, e.g. [16, 17]) and the correlation of their occurrence with the behavior of the material (contractive or dilative) and the test conditions (drained or undrained). Quite recently, the authors of the study mentioned in Ref [18] launched an experimental campaign in order to detect microscale precursors of visible laboratory-scale shear bands, observing diffusely distributed strain localization events, at the early stages of triaxial tests on sand, just few of which evolving toward visible shear bands.

In these studies, the most refined analyses have been obtained via X-ray computed microtomography, in particular high-speed synchrotron X-ray tomography, allowing the monitoring of non-steady fluid flow, on the one hand, and grain remodeling on the other. Only quite recently did the interest of the scientific community start to shed light on hydromechanical coupled localizations as preferential fluid paths in shear-banded regions in porous media, using neutron tomography techniques. In the study mentioned in Ref [19], it has been shown that an imbibition front accelerates when flowing through a compaction shear band. Conversely, the authors of the study mentioned in Ref [20] proved that the propagation of a non-wetting fluid front through a granular material can generate fracture, altering the local void ratio distribution: at the tip of the fracture, the grains are displaced and so the void ratio increases, while at the shoulder of the fracture, the grains are compacted and so a decrease in the void index is observed.

In this context, the present study aims at investigating the hydromechanical coupled instabilities in granular materials. The aim of the study is to understand the effect of a fluid front propagation on the skeleton remodeling, under biaxial loading, using a new biaxial setup endowed with a high-resolution optical system. Drainage experiments have been conducted by injecting air into a water-saturated sand sample. Then, the strains induced by air invasion are quantified using digital image correlation. The goal is not to conduct a full parametric study of the behavior of the sand assembly but to prove that the heterogeneous infiltration of gas through an initially saturated granular medium can possibly induce strain localization.

The article is organized as follows: Section 2 is devoted to describing the experimental setup; in Section 3, the results of the experimental campaign are illustrated, and in Section 4, a critical analysis of them is presented. In the last section, conclusions and perspectives of the study are briefly summarized.

2 Materials and Methods

To study the effect of hydraulic instabilities on the skeleton remodeling of a granular material, a new biaxial machine, adapted to partially saturated geomaterials and endowed with high-resolution cameras, has been designed at Ecole Centrale de Nantes. This is a suction-controlled biaxial equipment, derived from the classical Bishop–Wesley device, characterized by two couples of transparent windows, one directly confining the sample and the other one allowing to get visual access to the pressure chamber; the two couples of transparent windows are collocated on both sides of the equipment, which guarantees real-time visualization and recording of the occurring instabilities, using two cameras of 50 MPx, each equipped with telecentric lens (see Figure 1).


FIGURE 1. Biaxial machine endowed with its optical system.

2.1 Experimental Setup

The experimental setup allows applying both mechanical loading and hydraulic loading on a prismatic sample. The mechanical loading is given by a cell pressure provided by air compressors, which confines the sample cell inside the pressure chamber, and a vertical force/displacement provided by the thrust of the ram. The hydraulic loading allows developing suction-controlled tests or inducing a fluid flow through the sample, with a separate control of the fluid pressure at an injection and an extraction point placed at the bottom and top of the sample, respectively. In order to carry out tests in partially saturated conditions, a pressure–volume control for liquids and an air compressor for the gas to be injected are available. The acquisition and data control system is designed to manage a simultaneous multi-pressure control for three different fluids (gas or liquids) and the confining pressure. The maximum operating pressure is 1,600 kPa.

Therefore, thanks to this experimental setup, the fluid front penetration during a drainage or an imbibition process can be observed for different pairs of fluids (water–air, oil–air, water–oil, etc.). In this study, the focus is on drainage tests with air and water.

As previously mentioned, the sample cell of internal dimensions 40 × 50 × 11 mm3 differs from standard biaxial or triaxial cells because of two sapphire windows directly in contact with the specimen from the front- and backside views. In particular, the thickness of the sample has been chosen to be as small as possible to satisfy the conditions of two-dimensional fluid flow typical of the Hele-Shaw cell and to provide a minimal surface capable of sustaining the mechanical loading applied to the sample. The sample is not wrapped in any latex membrane. This geometry, mainly designed to enable the tracking of fluid propagation at the surface of the sample and to get full-field measurements of displacements, makes the sample’s sealing challenging. In order to face possible critical conditions, with respect to air intrusion into the sample, the lateral sides of the sample are put in contact with two latex membranes which directly transmit a slightly amplified cell pressure to them. This is obtained via a direct junction of the pressure chamber (labeled as (CO) in Figure 2) with a system made of burettes and a suitably defined pressure amplifier, which allows to apply a side pressure incremented of at least 10 kPa to the lateral side of the sample with respect to the cell pressure. From now on, all acronyms which we are referring to are mentioned in Figure 2. The top and bottom of the sample are confined by metallic blocks. They enable a uniform entry and exit of the in-place fluids during an upward or a downward flow through the porous stones embedded in them.


FIGURE 2. Biax block schematics provided by Megaris, Italy.

The optical system is constituted by a 50 MPx camera (VC-50MX) mounted with a telecentric lens of magnification 0.75. It covers a field of view 36.6 × 48.2 mm (6,004 × 7,904 pixels) at a resolution of 6.1  µm/pixel.

2.2 Procedure

The sample used is made of Fontainebleau sand NE 34 of a mean diameter d50 = 210 μm, which is considered as a reference material in many laboratory experiments and has been characterized in our laboratory in the study mentioned in Ref [21].

The drainage test is developed by injecting air into a sample previously saturated by deaerated water. In the following, the three preliminary steps of the drainage test are described in detail: the sample preparation phase, the confinement of the sample, and its saturation. The tests were conducted at room temperature.

2.2.1 Sample Preparation Phase

The preparation of the sample is a crucial phase to achieve a successful drainage test. The sample cell is constructed by assembling four parts: the top and bottom porous stones together with the left and right side membranes blocks. One side of the sample is then closed by the first sapphire window, then the sample block is placed horizontally to be filled with sand. Two removable confining brackets are put in place to fix the lateral bounds of the sample. The sample has been prepared by pluviation. However, the small dimensions of the sample, especially the thickness of 11 mm, make the sand-filling process challenging. In order to fit the geometrical constraints, due to the small dimensions of the sample, and the interaction between the sample and several delicate parts of the confining cell, a suitable pluviation setup has been developed. It provides a uniform and symmetrical rainfall of the sand and repeatable densities. Furthermore, this preparation method allows testing different densities by modifying either the falling height or the diameter of the holes of the sieve. Here, the apparent sand density is given and equal to 1.70 g/cm3 for a falling height equal to 250 mm; the hole diameter of the perforated plate is 2 mm; and the spacing between the holes 10 mm.

The backside is closed by the second sapphire window, and the specimen is kept for at least 24 h to dry the latex used at the corners of the sample to improve sealing efficiency. Then the sample cell is screwed to the chamber, and the scheme of the biaxial apparatus is given in Figure 2. Presaturation of all the hydraulic connections should be done, prior to the saturation phase, to reduce, as much as possible, the quantity of air present in the hydraulic circuit. The lateral left and right membranes are filled with silicon oil using a peristaltic pump (P). Once the membranes are saturated with oil, the brackets inside the membranes’ chambers, retaining the sample all over its height, are moved away and their confinement effect is replaced by the confining stress applied on the longitudinal sides of the sample via the cell pressure amplifier described in Section 2.1. The windows of the pressure chamber are closed, and the optical system is calibrated to film the center of the sample.

2.2.2 Confining Phase

An initial confining pressure of 45 kPa is applied to the drained sample, whereas the valve (FX), depicted in Figure 2, is opened to the atmospheric pressure. The applied cell pressure is delivered by the high-pressure input (HPIN) into the pressure chamber through the valve (CI for cell input), so providing the confinement of the sample. When the valve (CO for cell output) was opened, the cell pressure is transmitted to two burettes, the reference burette (RB) and the measurement burette (MB), to measure the volume changes of the sample. The (MB) burette is connected to the two longitudinal membranes via the side down block (SD) and transmits the confining pressure to the lateral membranes, after this has been incremented by at least 10 kPa by an external amplifier placed between SD and MB. The side pressure therefore acts on the lateral side of the sample with an overpressure with respect to that of the chamber pressure to guarantee a good sealing of the sample. The cell and side pressures are measured by two transducers CP and SP, respectively. When the cell pressure reaches its target value, a sealing check is done by turning the system to undrained conditions ((FX) is closed) and examining if any increase in pressure inside the sample occurs. If the sample does not show any leakage, the system is put back in drained conditions ((FX) is open).

2.2.3 Saturation Phase

The upward saturation of the medium starts through the circuit SI (for the saturation input) to SO (for the saturation output), connected to the pressure-controlled water reservoir. Due to a small pressure difference between the bottom and top of the sample (10 kPa), water fills the bottom porous stone and then percolates through the sample to saturate it. The water is then circulated through the sample to remove as much as possible the trapped bubbles of air in the circuit. To complete this saturation, the Skempton approach, detailed in the study mentioned in [22], is followed. A sequence of stepwise increments of the confining stress of Δσ = 50 kPa is applied (using the ramp mode of the machine controller) in undrained conditions, and the corresponding increment of the pore pressure ΔPp is measured via the SatP pressure transducer. At each stage of the loading path, the Skempton coefficient B, defined as the ratio between ΔPp and Δσ, is calculated considering it as an indicator of full saturation when it tends to 100%. Until the threshold is reached, once the pore pressure stabilizes, a back pressure is applied through the valve (FO) to maintain the same effective stress at the end of each stage. This valve is kept open until the injected water volume from the fluid reservoir toward the sample, measured by the sensor (FV), ceases. This process is repeated stepwise unless B becomes constant.

The different steps are summarized in Figure 3. The figure shows a sequence of cell pressure increases in blue, the resulting pore pressure measurement in red, the applied back pressure in yellow, and the consequent water injection in green. It is evident that the more steps are realized, the lower is the quantity of the injected water. This indicates the dissolution of air bubbles in water and consequently the approaching of full saturation.


FIGURE 3. Cell pressure, pore pressure, back pressure, and the difference between side and cell pressures on the left axis, and volume of water injected into the specimen on the right axis versus time, during Skempton phase.

2.2.4 Injection Phase

Finally, the injection phase consists of the invasion of air into the saturated sand under specific hydromechanical conditions. At the top of the sample, back pressure is maintained constant and equal to the value of the pore pressure Pp at the end of the saturation phase. At the bottom, the specimen is connected through the valves SX and GX to a gas block supplied by a high-pressure input reservoir. The gas in contact with the sample is initially at a pressure equal to the pore pressure. It uniformly occupies the bottom porous stone in such a way that when the injection is triggered, the gas source is not pointwise but distributed homogeneously across the sample cross section. This is an important point as the expected localization in the fluid flux should not be induced at the injection point but should manifest because of the loss of stability of the fluid front.

The gas pressure is then increased by constant mode to the target value Pg, to launch the gas propagation. In other words, a capillary pressure Pc = PgPp is imposed on the sample as the difference between the pressure of the non-wetting phase (the gas) and that of the wetting phase (water). The in-place water is expelled instantaneously toward the fluid reservoir, while the gas infiltrates the sample. The outflow volume is recorded via the sensor (FV). During the gas percolation, an image acquisition, at a frequency of 30 Hz, captures the drainage instabilities. In the following section, the images acquired from one front side of the sample are presented, and the obtained gas pattern is related to the outcomes of the biaxial acquisition system.

3 Results

3.1 Comparison of Two Drainage Tests

As already mentioned, the obtained results are far from constituting an exhaustive experimental campaign and just provide an example of the capabilities of the experimental apparatus to enlighten some interesting effects of microscale hydromechanical couplings on the response of the granular material at the scale of the laboratory sample. In the following paragraphs, the results of two test cases are presented whose characteristics are recorded in Table 1.


TABLE 1. Drainage tests under two different hydromechanical stress conditions.

In test (1), the imposed capillary pressure (30 kPa) is lower than the initial effective stress (40 kPa), provided by the difference between the total stress and the pore-water pressure in fully saturated conditions (Sr = 1). The gas invasion induced by the imposed capillary pressure implies a strongly ramified pattern. The fingered zone corresponding to the left side of the sample (see Figure 4), is the one which will be elaborated at a later point using digital image correlation.


FIGURE 4. Fontainebleau sand sample invaded by gas: on the left, drainage test (1); on the right, drainage test (2). The size of the image is height = 48.2 mm and width = 36.6 mm. The black rectangle marks the studied zone.

In test (2), the same imposed capillary pressure as in test (1) is considered, which, however in this case, is higher than the initial effective stress (15 kPa). Also, in this circumstance, a quite ramified pattern can be observed even if fingers are more easily identifiable all along the sample (see Figure 4).

The volume of the water expelled during gas injection is monitored by a volume sensor at the top of the sample, and the measure of the outflow can therefore be deduced. The results for the two tests are plotted in Figure 5: on the left, the volume change versus time, and on the right, the evolution of the outflow versus time. Assuming the product of air viscosity times the output flow, at time t = 1s, before the gas reaches the top of the sample, as a measure of the viscous forces, the logarithm of the capillary number is about log  Ca1 ≈ − 6.2 and log  Ca2 ≈ − 6.5 for test (1) and test (2), respectively, while the logarithm of the viscosity ratio is log  M = −2. A qualitative comparison with Lenormand’s [9] phase diagram implies that the two tests can be classified as in the intermediate zone between viscous and capillary fingering regions, test (2) being closer to the regime of capillary fingering than test (1). This is clearly a purely qualitative characterization of the flow regime as no account is taken of the presence of the granular skeleton.


FIGURE 5. Comparison of fluid volume expelled and outflow for the drainage tests (1) and (2).

We remark that the obtained patterns look qualitatively consistent with those of the study in Ref [23], test (2) exhibiting a regime more similar to capillary fracturing than test (1) which can be explained by the low confining load. However, no significant difference can be observed in the time the gas needed to attain the top of the sample in the two tests, which confirms that pore opening induced by gas injection into an initially saturated granular medium mainly depends on the release of capillary energy necessary to create a new fluid–fluid interface and only secondarily on the confining stress via porosity (or void index) current value (see [20]).

3.2 Image Processing and Strain Maps

In this section, a brief presentation of the method adopted to characterize the effect of fluid propagation on the remodeling of skeleton structure is provided. The main idea is to compare the gray level of the acquired images after gas invasion, called from now on the deformed state g, to the gray level of the reference image before gas invasion, labeled as the reference state f, using digital image correlation (DIC).

DIC is an image analysis technique that allows to explicitly calculate the displacement u of a pattern with respect to its reference state through time, based on the optical flow


for example, requiring that the gray level at a position x in the reference image should be found at a position x + u in the deformed image (the interested reader can find more details in [24]). In the present case study, the invasion of gas however modified the gray level of the medium between the reference water-saturated medium and the air-invaded partially saturated one, even without any deformation. To solve this issue, a correction of gray level is introduced, in the spirit of filtering out gas invasion-induced gray level changes.

The approach used for DIC is the so-called global approach, detailed in [25]. A finite element mesh of quadrilateral (square-shaped) elements is drawn over a region of interest (ROI). Eq. 1 is solved iteratively to get the nodal displacements u, taking into account a local and element-wise correction of gray level incorporated into the DIC algorithm that is called, zero-mean normalized sum of squared difference (ZNSSD), shown in Eq. 2:


where the index i = 1 denotes the deformed image and the index i = 0 denotes the reference image. At each iteration, the gray level of the reference image Im0 and the deformed image Im1 were corrected into Im0̃ and Im1̃ using the mean μi and standard deviation σi of the gray levels inside each reference element and its corresponding in the deformed state, respectively. This calculation is carried out using Ufreckles software [26].

In addition, a Tikhonov regularization with a cutoff wavelength of 15 pixels is used in order to act on the high-frequency displacement and to filter them out, so as to reduce the noise of the measurement [27].

Different mesh sizes have been tested to accurately capture the movements of grains. A mesh size of five pixels, which corresponds to 1/8 of the characteristic length of the grain speckle, has been finally retained as the most appropriate one, whereas coarser meshes, as for 10 or 20 pixels, are not capable of capturing grain movements, therefore hiding pore-scale deformation mechanisms.

The direct outcome of the algorithm is the nodal displacements; strains are then calculated element-wise via the element shape functions.

To analyze the impact of fluid invasion on the solid matrix, the first invariant of the strain tensor, the volumetric strain ϵv, is computed as follows:


where ϵxx is the horizontal strain and ϵyy is the vertical strain. ϵv indicates the volume change of the material: when positive, it describes swelling, and when negative, it describes shrinkage. The volumetric strain maps for a series of deformed images of drainage tests (1) and (2) are shown in Figure 6.


FIGURE 6. Volumetric strain maps for a mesh size of five pixels and a Tikhonov regularization of 15 pixels corresponding to a series of images captured during tests (1) and (2). These time shots are marked by * in Figure 9 and Figure 10.

The error on the volumetric strains due to the sensitivity of the sensors to vibrations and possibly induced by the change of contrast due to the gray-level correction procedure has been quantified by testing the DIC algorithm on a set of fixed experimental images and a set of images virtually saturated with a gray-level distribution extracted from that of the fully saturated sample, respectively. The order of magnitude of strains referred to the adopted mesh size is in the two cases of 1% and 0.9%.

The black line which identifies the interface between the wetting and the non-wetting phase in the considered sequence of maps has been identified using a phase-field approach inspired by a gradient damage model. At the base, this method aims to replace the discontinuity of a medium such as a fracture, with a continuous field d spread over a distance lc that describes the smeared crack ([28]; [29]; [30]). Thinking of the application to partially saturated images, the part of the images invaded by air was considered as the ≪ fractured ≫ zone, corresponding to d different from zero, while the remaining part is considered undamaged (d = 0). The phase field d was therefore inversely proportional to the degree of saturation of the medium. The lc parameter, in this case, described the thickness of the partially saturated transition zone. For a level set value d = 0.25, contour lines can be plotted to separate the drained zone from the undrained one. The obtained interfaces for the sequence of images of tests (1) and (2) are presented in Figure 7.


FIGURE 7. Series of raw images acquired during the gas invasion for test (1) in the first row and test (2) in the second row. The timestep between the deformed images is 1/3s for test (1) and 2/3s for test (2). The left column represents the reference states for both tests, and the remaining columns are for the deformed ones. The bright color indicates the gas, while the dark color indicates the water. The black lines represent the interfaces between saturated and unsaturated part of the sample.

4 Discussion

Following the drainage process via the optical system described in Section 2.1 allows tracking the fluid front and getting the full-field strain maps as described in Section 3.2. In this section, a comparison of the results of the two tests described in Section 3 is discussed.

As already mentioned, both the considered test cases are associated to quite ramified patterns that are formed by gas injection at the bottom and quickly propagate toward the top of the sample. The analyses discussed in what follows concern just the two subdomains, identified with the two black boxes of Figure 4, of the complex patterns which show up when the imposed capillary pressure is smaller or greater than the initial effective stress, respectively.

In both tests, the evolution of the fluid front is accompanied by the simultaneous appearance of narrow clusters of elements associated with positive and negative volumetric strains, in the drained zone of the granular medium. This can be observed assessing a close-up window at the tip of the finger between two consecutive instants, as shown in Figure 8. When the interface progresses between ti and ti + 0.033s, two clusters of positive and negative strains show up simultaneously; the order of the magnitude of the measured strains is clearly larger than the previously quantified errors.


FIGURE 8. Close-up windows on volumetric strain maps at the tip of the finger at ti and ti + 0, 033s during the test (1), obtained using a mesh size of five pixels and a Tikhonov regularization of 15 pixels.

Considering the characteristic size of the finite elements used to develop the DIC analysis, these strain patterns could be interpreted as a straightforward consequence of grain displacements: in this spirit, the local shrinkage corresponds to grain entanglement, for instance, grain reorganization from an initial configuration toward a denser one and consequently negative volumetric strain; local swelling to grain disentanglement, for instance, grain extrication giving rise to a looser configuration and consequently positive volumetric strain. No significant induced strain has been observed in the back of the non-wetting/wetting interface, when the complete drainage of the pore network has been attained: strains once triggered remain (almost) constant in time during fingering propagation, as shown in Figure 6.

A speculative but illustrative explanation of the relation between the volumetric strain distribution, observed in Figure 8, and the saturation evolution at the grain scale during drainage, resides in the pore-scale mechanisms that characterize the drainage phenomenon, and specifically the Haines jump events. A Haines jump event consists of a sharp drop of the capillary pressure when the gas enters the target pore which is accompanied by the imbibition of the neighbor channels. During drainage, the non-wetting phase infiltrates first the largest channel as, according to the Laplace law, it has the lowest entry pressure. When the meniscus passes from the pore neck into the pore, initially occupied by the wetting phase, the capillary pressure drops down because of the change in meniscus’ curvature.

Considering the classical definition of the Bishop effective stress, valid in the regime of partial saturation, for instance, σ′ = σ + PcSr, the local effective stress drops down too, which explains the local swelling of the granular medium (positive volumetric strains). In the adjacent channels, however, where the wetting phase pushed the non-wetting fluid to sustain the occurring drainage process, the corresponding rapid imbibition is accompanied by an increase in the capillary pressure, as explained in the study mentioned in Ref [31]; as a result, an increase in the local effective stress occurs which is interpreted as a shrinkage of the porous network (negative volumetric strains).

Even if neither a local measure nor a fluctuation of the global capillary pressure can be detected via the experimental setup, the qualitative correlation between local strain distribution and gas front propagation seems to corroborate such an interpretation. Sun and Santamarina [10] found similar results for the deformation of capillary ducts during drainage.

The time evolution of the average strains within the finger formations identified in Figure 4 has also been evaluated for the two test cases. In Figure 9 and Figure 10, the average positive and the absolute value of the average negative volumetric strains, panel (A); the standard deviation of positive and negative volumetric strains, panel (B); the frequency of positive and negative volumetric strains within the finger, panel (C); and the average volumetric strain within the finger, panel (D), are depicted. It is evident that as the average positive and the absolute values of the average negative volumetric strains are almost the same all along the fingering nucleation and propagation process in both tests, this is not the case for the average volumetric strain which is positive in both cases but higher and with a more negative slope in test case (1) than in test case (2). This result is corroborated by the frequency of positive volumetric strains, within the finger, higher than that of negative ones.


FIGURE 9. Time evolution of mean, standard deviation, and frequency of positive and negative volumetric strains captured inside the finger (A, B) and (C), respectively, and mean of volumetric strains inside the finger (D) for test (1). The time instants corresponding to the images of test (1) in Figure 7 are labeled by the symbol *. The reference t0 corresponds to the instant of gas entering the sample.


FIGURE 10. Time evolution of mean, standard deviation, and frequency of positive and negative volumetric strains captured inside the finger (A, B) and (C), respectively, and mean of volumetric strains inside the finger (D) for test (2). The time instants corresponding to the images of test (2) in Figure 7 are labeled by the symbol *. The reference t0 corresponds to the instant of gas entering the sample.

In other words, the higher the initial effective stress, and therefore the confinement, is, the more relevant the initial positive strain associated to the fluid–fluid interface nucleation will be. Moreover, as the limit value of the volumetric strain at the stationary state looks similar in tests (1) and (2), the higher the effective stress is, the more important the strain decay will be.

These results indeed underline the double-scale behavior of strains, induced by drainage, in granular media: an intrinsic “discrete nature” of strains at the grain scale, associated with grain displacements, and a “continuous behavior” of strains at the mesoscopic scale of the sample, where a non-vanishing cumulated strain is appreciated as a result of its spatial distribution at the grain scale.

5 Conclusion

In this study, the first results obtained by exploiting the new biaxial apparatus available at Ecole Centrale de Nantes for the characterization of micro-scale hydromechanical couplings in granular media have been presented. In particular, evidence of grain remodeling induced by drainage in the regime of capillary fingering has been found, and the corresponding strains have been quantified via DIC techniques, once the front of the invading draining fluid is identified. To our knowledge, this is the first time where this coupling is observed and quantified.

Results presented in the article do not aim at providing a full parametric study of the considered sand assembly but to give a proof-of-concept of the main idea that heterogeneous infiltration of gas through an initially saturated granular medium can be responsible for local strain concentration, which indeed has been observed in the two tests presented in Section 3 and discussed in Section 4.

Further developments will be devoted to validate the proposed interpretation of the relation between grain remodeling and Haines jump events, driving the gas injection via an imposed flow rather than a given gas pressure; this will allow for detecting possible fluctuations of the capillary pressure in a similar way as observed by [12]. Moreover, with the aim of considering possible applications of these results to the characterization of the tightness of the sealing caprock of aquifer rock reservoirs, a similar analysis will be developed considering materials of much smaller intrinsic permeability or hydraulic conductivity. In this case, the use of fluids other than water and air will be recommended.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

RA, GS, and JR contributed to conception and design of the study. RA performed the tests and analyzed the results. All authors contributed to manuscript revision, read, and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.


The authors would like to acknowledge the support of the French National Research Agency (ANR), project STOWENG (Project-ANR-18-CE05-0033). The authors acknowledge Megaris company (Adolfo Cavallari and his team) for the realization of Biax the new biaxial machine and its installation at Centrale Nantes.


WP, wetting phase; NWP, non-wetting phase; Δσ, increment of confining pressure (kPa); ΔPp, increment of pore pressure (kPa); B, Skempton coefficient (%); Pp, pore pressure (kPa); Pg, gas pressure (kPa); Pc, capillary pressure (kPa); Sr, degree of saturation; σ, confining pressure (kPa); σ0, initial effective stress (kPa); σ′, effective stress for partially saturated conditions (kPa); f, gray level in the reference state; g, gray level in the deformed state; u, nodal displacement vector; Imi,, gray level for reference image if i = 0, for deformed images if i = 1; μi, mean value of gray levels; σi, standard deviation of gray levels; d, phase-field parameter.


1. Iglauer S, Paluszny A, Pentland CH, Blunt MJ. Residual Co2 Imaged with X-ray Micro-tomography. Geophys Res Lett (2011) 38, L21403. doi:10.1029/2011gl049680

CrossRef Full Text | Google Scholar

2. England WA, Mackenzie AS, Mann DM, Quigley TM. The Movement and Entrapment of Petroleum Fluids in the Subsurface. J Geol Soc (1987) 144:327–47. doi:10.1144/gsjgs.144.2.0327

CrossRef Full Text | Google Scholar

3. Xu L, Davies S, Schofield AB, Weitz DA. Dynamics of Drying in 3d Porous media. Phys Rev Lett (2008) 101:094502. doi:10.1103/PhysRevLett.101.094502

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Buijze L, Guo Y, Niemeijer AR, Ma S, Spiers CJ. Nucleation of Stick‐Slip Instability within a Large‐Scale Experimental Fault: Effects of Stress Heterogeneities Due to Loading and Gouge Layer Compaction. JGR Solid Earth (2020) 125:e2019JB018429. doi:10.1029/2019JB018429

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Verberne BA, van den Ende MPA, Chen J, Niemeijer AR, Spiers CJ. The Physics of Fault Friction: Insights from Experiments on Simulated Gouges at Low Shearing Velocities. Solid Earth (2020) 11:2075–95. doi:10.5194/se-11-2075-2020

CrossRef Full Text | Google Scholar

6. Hall SA, Bornert M, Desrues J, Pannier Y, Lenoir N, Viggiani G, et al. Discrete and Continuum Analysis of Localised Deformation in Sand Using X-ray μCT and Volumetric Digital Image Correlation. Géotechnique (2010) 60:315–22. doi:10.1680/geot.2010.60.5.315

CrossRef Full Text | Google Scholar

7. Spetz A, Denzer R, Tudisco E, Dahlblom O. Phase-field Fracture Modelling of Crack Nucleation and Propagation in Porous Rock. Int J Fract (2020) 224:31–46. doi:10.1007/s10704-020-00444-4

CrossRef Full Text | Google Scholar

8. Spetz A, Denzer R, Tudisco E, Dahlblom O. A Modified Phase-Field Fracture Model for Simulation of Mixed Mode Brittle Fractures and Compressive Cracks in Porous Rock. Rock Mech Rock Eng (2021) 54:5375–88. doi:10.1007/s00603-021-02627-4

CrossRef Full Text | Google Scholar

9. Lenormand R. Liquids in Porous media. J Phys Condens Matter (1990) 2:SA79–SA88. doi:10.1088/0953-8984/2/s/008

CrossRef Full Text | Google Scholar

10. Sun Z, Santamarina JC. Haines Jumps: Pore Scale Mechanisms. Phys Rev E (2019) 100:023115. doi:10.1103/PhysRevE.100.023115

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Eriksen FK, Toussaint R, Måløy KJ, Flekkøy EG. Invasion Patterns during Two-phase Flow in Deformable Porous media. Front Phys (2015) 3:81. doi:10.3389/fphy.2015.00081

CrossRef Full Text | Google Scholar

12. Berg S, Ott H, Klapp SA, Schwing A, Neiteler R, Brussee N, et al. Real-time 3d Imaging of haines Jumps in Porous media Flow. Proc Natl Acad Sci (2013) 110:3755–9. doi:10.1073/pnas.1221373110

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Bultreys T, Boone MA, Boone MN, De Schryver T, Masschaele B, Van Loo D, et al. Real‐time Visualization of H Aines Jumps in sandstone with Laboratory‐based Microcomputed Tomography. Water Resour Res (2015) 51:8668–76. doi:10.1002/2015wr017502

CrossRef Full Text | Google Scholar

14. Zacharoudiou I, Boek ES. Capillary Filling and haines Jump Dynamics Using Free Energy Lattice Boltzmann Simulations. Adv Water Resour (2016) 92:43–56. doi:10.1016/j.advwatres.2016.03.013

CrossRef Full Text | Google Scholar

15. Zacharoudiou I, Boek ES, Crawshaw J. Pore-scale Modeling of Drainage Displacement Patterns in Association with Geological Sequestration of Co 2. Water Resour Res (2020) 56:e2019WR026332. doi:10.1029/2019wr026332

CrossRef Full Text | Google Scholar

16. Desrues J, Chambon R, Mokni M, Mazerolle F. Void Ratio Evolution inside Shear Bands in Triaxial Sand Specimens Studied by Computed Tomography. Géotechnique (1996) 46:529–46. doi:10.1680/geot.1996.46.3.529

CrossRef Full Text | Google Scholar

17. Lenoir N, Bornert M, Desrues J, Bésuelle P, Viggiani G. Volumetric Digital Image Correlation Applied to X-ray Microtomography Images from Triaxial Compression Tests on Argillaceous Rock. Strain (2007) 43:193–205. doi:10.1111/j.1475-1305.2007.00348.x

CrossRef Full Text | Google Scholar

18. Desrues J, Andò E, Bésuelle P, Viggiani G, Debove L, Charrier P, et al. Localisation Precursors in Geomaterials? In: International Workshop on Bifurcation and Degradation in Geomaterials. Cham: Springer (2017). p. 3–10. doi:10.1007/978-3-319-56397-8_1

CrossRef Full Text | Google Scholar

19. Hall SA. Characterization of Fluid Flow in a Shear Band in Porous Rock Using Neutron Radiography. Geophys Res Lett (2013) 40:2613–8. doi:10.1002/grl.50528

CrossRef Full Text | Google Scholar

20. Shin H, Santamarina JC. Fluid-driven Fractures in Uncemented Sediments: Underlying Particle-Level Processes. Earth Planet Sci Lett (2010) 299:180–9. doi:10.1016/j.epsl.2010.08.033

CrossRef Full Text | Google Scholar

21. Yin K, Fauchille AL, Di Filippo E, Othmani K, Branchu S, Sciarra G, et al. The Influence of Mixing Orders on the Microstructure of Artificially Prepared Sand-clay Mixtures. Adv Mater Sci Eng (2021) 2021. doi:10.1155/2021/8552224

CrossRef Full Text | Google Scholar

22. Head KH. Manual of Soil Laboratory Testing, Vol. 3. Chichester, UK: John Wiley & Sons (1998). p. 2. effective stress tests.

Google Scholar

23. Holtzman R, Szulczewski ML, Juanes R. Capillary Fracturing in Granular media. Phys Rev Lett (2012) 108:264504. doi:10.1103/physrevlett.108.264504

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Sutton MA, Orteu JJ, Schreier H. Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications. Boston, MA: Springer Science & Business Media (2009).

Google Scholar

25. Besnard G, Hild F, Roux S. "Finite-Element" Displacement Fields Analysis from Digital Images: Application to Portevin-Le Châtelier Bands. Exp Mech (2006) 46:789–803. doi:10.1007/s11340-006-9824-8

CrossRef Full Text | Google Scholar

26. Réthoré J. Ufreckles (2018). URL 1433776.

Google Scholar

27. Witz J-F, Réthoré J, Hosdez J. Regularization Techniques for Finite Element Dic. In: International Digital Imaging Correlation Society. Cham: Springer (2017). p. 137–40. doi:10.1007/978-3-319-51439-0_33

CrossRef Full Text | Google Scholar

28. Bourdin B, Francfort GA, Marigo J-J. The Variational Approach to Fracture. J Elasticity (2008) 91:5–148. doi:10.1007/s10659-007-9107-3

CrossRef Full Text | Google Scholar

29. Miehe C, Hofacker M, Welschinger F. A Phase Field Model for Rate-independent Crack Propagation: Robust Algorithmic Implementation Based on Operator Splits. Comput Methods Appl Mech Eng (2010) 199:2765–78. doi:10.1016/j.cma.2010.04.011

CrossRef Full Text | Google Scholar

30. Nguyen TT. Modeling of Complex Microcracking in Cement Based Materials by Combining Numerical Simulations Based on a Phase-Field Method and Experimental 3D Imaging. Paris Est: Université Paris-Est (2015). Ph.D. thesis.

31. Lenormand R, Zarcone C, Sarr A. Mechanisms of the Displacement of One Fluid by Another in a Network of Capillary Ducts. J Fluid Mech (1983) 135:337–53. doi:10.1017/s0022112083003110

CrossRef Full Text | Google Scholar

Keywords: drainage, hydraulic instability, biaxial apparatus, fingering, solid remodeling, partially saturated medium

Citation: Al Nemer R, Sciarra G and Réthoré J (2022) Drainage Instabilities in Granular Materials: A New Biaxial Apparatus for Fluid Fingering and Solid Remodeling Detection. Front. Phys. 10:854268. doi: 10.3389/fphy.2022.854268

Received: 13 January 2022; Accepted: 24 February 2022;
Published: 18 March 2022.

Edited by:

Ferenc Kun, University of Debrecen, Hungary

Reviewed by:

Hiroaki Katsuragi, Nagoya University, Japan
Bjornar Sandnes, Swansea University, United Kingdom

Copyright © 2022 Al Nemer, Sciarra and Réthoré. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Rana Al Nemer,