*

e experience of pursuing a PhD could have been expressed using several words. In view of being concise, here I choose one: ful lling.

ere are several individuals that have accompanied me in this journey of over three years. While a single typeset page would not do justice to their contributions, I owe them all a sincere acknowledgment.

Foremost, I would not have been able to complete this thesis without the unwavering guidance of my advisers: Giulio Sciarra and Panagiotis Kotronis. I express my deepest gratitude to both of you for your continued support and excellent direction during my rst steps. Giulio, thank you for introducing me to this strangely complex world of porous media and guiding me onto the right path whenever I tend to meander. I especially admire and envy your unrelenting patience that has been constructive in shaping my thoughts and work. Panos, I am indebted to you for your enthusiasm and your meticulous comments on my work. Many a time, your encouragement has pushed me to pursue aspects that I originally deemed arduous. I am grateful for the jury members and deeply honored by their kind acceptance of the invitation. I thank M. Sergey Gavrilyuk and M. Pierre Seppecher for accepting to spend their invaluable time in reading my thesis. I thank also the commi ee that has followed-up my work for this thesis: M. Corrado Maurini and M. Pierre Bésuelle, whose guidance has been indispensable.

I have been fortunate to be part of the GeM laboratory especially because of the environment it has created to nourish younger researchers at the beginning of their career. I thank M. omas Heuzé and M. Laurent Stainier for giving me the opportunity to teach along side them which has been a rewarding experience. I thank all my colleagues for being no less than kind towards me and for providing plentiful positive energy. I am grateful to Sara and entin who have shared the o ce with me all these years. Sara, I thank you deeply for the friendship and warmth that you have shared throughout.

entin, I thank you for sharing your enthusiasm and, mostly french, anecdotes. I thank, Katia for her energetic presence and support, Eric for sharing his know-how and for keeping the Titan cluster running through the heaviest of loads. e colleagues in the ANR-STOWENG project deserve a special mention: Rana for her meticulous experimentation, Mohammed for sharing the other side of modeling and colleagues of the 3SR laboratory at Université Grenoble Alpes for sharing their experimental observations. Special thanks goes to people who gathered around the table for lunch everyday sharing their laughs: Alex, Georges, Praveen, Mickäel, Filippo, Diego, Jinane, Imane/s, Youmna, Ahmad, Abdallah, Guilleme e, Aliénor, … the list is long but is incomplete without mentioning Julie a for bringing much joy. I thank my dear friend Vasu for always being an example and bearing with my peculiar humor throughout our long journey together. A heartfelt thanks goes to Serge Lefrant, Annie Simon and to Cataldo, Clara for being kind enough to open their doors to me and to Claude… for being Claude! I am indebted to the numerous individuals who have spent their hours in supporting various open-source and community-driven so wares/libraries that have been invaluable for my thesis: special mention goes to the Linux-Mint team, FEniCS, PETSc/TAO and SLEPc.

I nally thank my large family, the list of which would ll another page. My mother and grand mother, for all their sacri ces and patience throughout my endeavors. My grandfather, for encouraging me to explore the world and whose presence I dearly miss. Both my uncles and my aunts who have been the pillars of my growth and my siblings who give me a sense of togetherness. is page would be incomplete without thanking my most beloved one, Graziana, for being an unrelenting support and source of energy for me all these years. Your essential presence has been THE most exceptional in uence through this journey.

All these individuals and many others, known and unknown, have had their supportive hand in my journey up until now and hopefully in the future.

*

Hydraulic and mechanical instabilities in geomaterials refer to a variety of non-linear phenomena that can be triggered by heterogeneities inherent to such materials. Hydraulic instabilities in partially saturated conditions can manifest themselves as heterogeneous uid invasion causing ' ngering' phenomenon. Mechanical instabilities on the other hand can present themselves as strain localizations and/or fractures. ese instabilities and their associated coupling pose a major obstacle for a myriad of applications involving geomaterials like Carbon dioxide (CO 2 ) sequestration, rain water in ltration and also for reliable predictions such as for contaminant ow in ground waters. Existing classical models do not resolve this behavior due to their lack of pa ern-forming ingredients in their formulation and thus being stable against perturbations. e essence of current thesis work is to propose and investigate modeling techniques that allow to describe these instabilities in a simple and natural manner. e constitutive approach adopted in this thesis is that of micro-structured continua, in particular that of enhanced continua with a constitutive law depending on the gradient of so-called phase eld variables.

e rst part of the work deals with description of a uid-uid front that has been modeled as a di used interface. is has been done by interpreting the presence of two uids within the pore space as a single non-uniform uid and the degree of saturation of one of the uids as the corresponding phase eld. While classical one-to-one relation between capillary pressure and saturation degree is expected to describe the retention properties of the porous network due to its texture, an enhanced relation is obtained by prescribing a chemical potential in the spirit of Cahn-Hilliard type modeling of multi-phase uids. is enhanced relation together with a non-local energy contribution provides the ingredients required to describe uid-ngering and non-local phenomena such as pinching and coalescence of uids due to variations in capillary forces. Using asymptotic matching techniques it has been shown analytically, in a one-dimensional se ing, that the particular choice of energy contributions employed in the above framework allows to replicate 'overshoot' type non-monotonic saturation pro les during constantrate in ltration into relatively dry medium. is has been found to match qualitatively with experimental observation. Further, a slight non-convexity introduced into the ux function has been shown numerically to allow the modeling of drainage fronts, besides imbibition, without employing any additional complexities into the model. A careful linear stability analysis (LSA) of the homogeneous solutions against arbitrary perturbations has been done. e growth in time of the amplitude of imposed transverse perturbations is understood as the ngering phenomenon. e predictions of the LSA have been followed-up by two-dimensional simulations portraying the ability of the proposed model to describe uid ngering and segregation.

In the second part the triggering of a fracture within a drying porous medium has been studied. A prevailing modeling perspective, in the spirit of gradient damage modeling, has been rst tested for its ability to replicate periodic fracture formation as observed in representative experiments. Further an investigation is done, starting from the gradient damage modeling framework, but interpreting the presence of a fracture within a drying porous material as a loss of its capillary properties, thus allowing the passage of non-we ing uid under vanishing capillary pressure. is paradigm is v of particular interest in modeling cohesion-less and unconsolidated negrained soils, where the resistance against tensile loading is negligible and thus fracturing induced due to development of tensile stresses is not the prevailing phenomenon. Starting from the principles of variational approach used to construct the model, it has been shown that for su ciently strong desiccation, damage initiates homogeneously on the drying face while progressing into the body with time. e possible occurrence of bifurcations of this base solution has been analyzed, again in the framework of LSA.

is works sets the stage for various possibilities, the most natural one being the study of coupling between the above mentioned instabilities. Apart from that, two-dimensional simulations of uid-ngering have shown that the model predicts additional features of unstable ow, such as pinching and coalescence of the we ing phase, which need to be investigated using carefully designed experiments. Initiation of damage induced due to evolving drainage nger is also of particular interest in the context of earlier mentioned applications. Lastly, advanced numerical techniques can be sought a er for resolution of the above problems with an intent to provide accurate solutions more e ciently. It comes as a no surprise that since the industrial revolution human activity is a ecting, at a grand scale, the natural cycles that preserve the habitability of our planet. For instance, anthropogenic contribution of carbon dioxide (CO 2 ), especially due to the use of fossil fuels, to the carbon cycle is a major cause of its accumulation in the Earth's atmosphere, dropping pH of the oceans and increasing global temperatures. Each of these in turn have an array of devastating a ects on the delicate balance that is needed . Our accumulated knowledge is su cient to know that even though there exist sinks within this cycle that can potentially reduce the harmful impact, there is no existing natural savior that can act as a silver bullet (F ., ). However, various techniques have been proposed, a combination of which could do this job. Carbon Capture & Storage (CCS) being one of them, has proven (L B , ) over the past 50 years to be a promising solution especially because it ensures a smooth transition to renewable resources instead of a sharp shi in the energy supply. It is especially a ractive because it can be employed to directly target CO 2 generated at coal or gas red power plants and cement production plants, which account for over 40% of the worldwide emissions. is is because these plants act as stationary and predictable sources of CO 2 . Apart from CCS, a solution that focuses on the convergence towards a low carbon economy and balancing of supply and demand is the large-scale underground storage of energy. is has shown (P C S , ; M ., ) the potential to integrate renewable energy sources while e ciently managing the grid, thus increasing energy security.

While the operational aspects of the above mentioned solutions are similar, for the purpose of being concise we contextualize CCS. is involves capturing at a large scale the atmospheric CO 2 and then sequestering it at its super-critical conditions. At these conditions the gas retains its low viscosity but a ains an increased density, reducing substantially the required storage volume. At depths of around 1km and higher within the Earth's crust the conditions are just right for this purpose, thus providing a relatively vast storage volume. e injection of compressed CO 2 in high permeable rock formations at suitable geological sites involves displacement of the resident we ing uid, usually brine water (in saline aquifers) or residual oil (in depleted oil elds). A low permeable cap-rock usually acts as a natural barrier to trap the injected CO 2 , Fig. ( . ). As any other operation involving geomaterials, the above process comes with a myriad of complexities and risks. One of the major concerns is the proper containment of CO 2 during the injection and ensuing migration of the gas plume. Geomaterials are known to be inherently inhomogeneous with existing fractures and new ones on the cusp of being triggered. It is known and well-studied (W Z , ; C

., ; E S , ) that the previously dormant faults within cap-rocks can be triggered due to pressure changes within the vicinity. Another concern is the predictability of the ow itself. Any instability relative to an advancing uid front could lead to reduced e ciency of the operation (Y ., ; Z .,

). It can even pose an environmental risk if preferential paths for the uid ow are triggered within the cap-rock due to fault planes, thus providing an escape route to the rising or migrating CO 2 plumes (M ., ; G ., ; S Z , ). e motion of two immiscible uids, in the above case CO 2 and brine, one displacing the other through an inhomogeneous material which itself could fail is a coupled problem. While such coupled phenomena are contextualized above for CCS operations, there are numerous other applications within soil hydrology that fall victim to inhomogeneous uid ow and coupled failure mechanisms. Transport of solutes and contaminants to the ground waters, Fig. ( . ), is a phenomenon that has far-reaching consequences (G ., ; J , ; C ., ). While natural rain and irrigation waters replenish the ground waters, they also tend to carry the surface pollutants such as fertilizers and toxic wastes. However, the soil layers act as natural lters of these pollutants before they could reach the water table .  e e ciency of this process and the ability to make predictions of ground water contamination are known to be negatively impacted by the formation of preferential pathways due to unstable in ltration.

Another context where immiscible uid ow can trigger pore network failure is in the hydrate stability zone (HSZ) found under ocean sediments along continental margins (S , ). In HSZs it has been observed that methane hydrates, methane gas and brine can co-exist (G ., ). While methane hydrates are widely believed to be a potential source of energy supply, the prevailing risk is the release of trapped methane gas through existing or newly formed faults within the sediments. It has been hypothesized that continued gas accumulation below the HSZ can result in activation of fractures (Z S , ; L F , ) which could also be the reason for underwater oceanic land slides (P ., ) in addition to the release of the highly potent green-house gas.

M

While the prediction and control capabilities in real-world applications of the immiscible uid-uid displacement coupled to the porous matrix could be ridden with uncertainties, e cient modeling of the phenomena involved gives a head-start. Models provide new insights into individual phenomena involved within the coupling and their associated failure mechanisms. is has been the rationale for an ever-growing body of research within which the current thesis is a small part.

O

Failure in the current context is the formation of preferential pathways for the uids to ow. For immiscible uid-uid displacements this could be viewed as an instability of the moving front also called uid ngering. On the other hand, mechanical failure of the porous matrix that creates a pathway could be viewed as a localized re-structuring of the soil grains which could lead to an opening mode fracture or at least a localized modi cation of the material properties such that one uid can displace the other with ease. While classical theories that have been used over the years to model bi-phasic ow in deformable porous media are able to model ideal/stable processes, they have been known to be lacking the ability to model unstable phenomena and coupling involved. However, various modeling extensions have been proposed to remedy this lacking, each building upon the existing classical models. It is the purpose of the current thesis to propose such extended models and investigate their ability to describe both hydraulic and mechanical instabilities.

P

e current thesis is organized as follows. Chp.( ) is dedicated to the phenomenology, recounting the existing experimental investigations of the unstable phenomena. Particular observations are done of the evolution structure especially of the spatial uid distribution, that seem to trigger the said instabilities. Chp.( ) sets the stage for modeling with a concise recall of the classical model of partial saturation and a phase eld approach. Part II is dedicated to the modeling of hydraulic instabilities starting from the governing equations of the phase eld approach in the context of nondeformable porous network. An in-depth analysis of the one-dimensional solution structure is presented followed by a linear stability analysis under imposed perturbations and two-dimensional simulations. Part III concerns the mechanical instabilities within a deformable porous network. First an application of an existing technique to model fracture initiation is extended to the context of classical partial saturation. Following this, an investigation of an alternative triggering phenomenon of opening mode fractures, that is motivated by experimental observations, is investigated by revisiting the driving force employed in the modeling.

A

e current thesis study was conducted under the framework of a doctoral contract supported by the French ministry: Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation through the doctoral school: École doctorale Bretagne Loire -Science pour l'ingenieur. Also the reaserch work falls under the broad framework of the ANR-STOWENG project (ANR--CE -) (ANR STOWENG,

). e work in its entirety was conducted at the Institut de Recherche en Génie Civil et Mécanique (GeM), UMR CNRS, École Centrale de Nantes. Unstable uid ow and fracturing are two distinct phenomena that characterize failure of porous media. In this chapter we present evidence of these phenomena from the existing literature. e various regimes of ow morphology and the analogous nature of bi-phasic porous media ow to that of Hele-Shaw ow are highlighted. e role that gravity plays in in ltration and drainage scenarios is as well discussed with a focus on saturation pro les and their link to destabilization of the ow. e phenomenon of fracturing is presented with a prevailing treatment for opening mode fractures during drainage scenarios such as desiccation. An alternate phenomenological treatment is recounted, supported by experimental evidence. Overall, this exposition sets the stage for phenomenological motivations and qualitative evaluations of the modeling results.

. Stability of a moving interface between two immiscible uids of contrasting viscosities and densities is a classical problem in uid mechanics. Further when a solid substrate is involved the problem becomes more complex with the addition of we ability of the solid surface with respect to the uid combination. is is almost always the case when two uids are involved such as in the Hele-Shaw ow and in multi-phase ow in porous media such as soils. ) depicting qualitatively the zones of stable and unstable porous media ow. Irrespective of the we ing properties of the invading and defending uids these extremes of ow morphology were demonstrated experimentally to exist at similar extremes both for imbibition and drainage.

When the resident/ defending uid preferentially wets the solid substrate the displacement of it by a non-we ing uid is termed as drainage. e contrary is termed imbibition, when the invading uid is we ing. In both these scenarios, numerous experimental campaigns have been carried out by various researchers with an intent to characterize the unstable displacement and the ensuing pa ern formation. One can refer to the seminal works of Lenormand and his colleagues in the s (L , ; L ., ; L , ) and the other researchers who followed this path. See C . ( ); G A ( ), to name a few.

eir work has laid the foundations for a more systematic way to identify various regimes of such ow by employing the so-called 'Phase diagram' (L , ; L ., ). See Fig. ( . ). In essence, a stable and two morphologically di erent unstable regimes have been characterized by the dominant forces at play within the ow. ese forces are either the capillary forces or the viscous forces. ese forces were quanti ed by two non-dimensional numbers: capillary number, Ca, and viscosity ratio, M, de ned as,

Ca = V i µ i γ , M = µ i µ d , ( .
)
where γ is the surface tension between the two uids, µ and V are respectively the dynamic viscosity and mean velocity of the uids. Subscripts i and d represent the invading and defending uids respectively. e phase diagrams were plo ed with the log-scale of these two numbers. In essence the capillary number quanti es the ratio of viscous forces within the invading uid to that of the capillary forces acting at the various uid menisci. e viscosity ratio on the other hand quanti es the ratio of viscous forces within the invading and defending uids.

As summarized in L ( ) the stable regime is the one in which the viscous forces within the invading uid dominate both the viscous forces within the defending uid and the capillary forces. Viscous ngering regime is characterized by dominating viscous forces within the defending uid compared to those within the invading uid and the capillary forces. In this regime at initial times the front destabilizes into several ngers followed by some ngers growing faster than others and inhibition of shorter ngers. All the ngers tend to orient towards the exit forming a characteristic tree-like pa ern (L ., ). Capillary ngering regime has capillary forces dominating those due to viscous e ects and is characterized by ngers that grow in all the directions and even forming loops that trap the displaced uid. We ing conditions have been known to impact how the invading uid advances at the pore-scale and even for the same we ing condition, the pore-scale geometric features, such as the pore-to-throat size ratio in the case of imbibition and throat size distribution in the case of drainage, have consequential e ects on the invasion morphology. Notwithstanding these di erences, the phase diagrams revealed that these regimes of stable and unstable ow do exist at similar extremes in both imbibition and drainage.

. . E ect of gravity

It is to be noted that in the works of Lenormand and others cited above, the ow is within a slow regime (no inertial e ects) and in the absence of gravity. Now in a practical se ing, gravity almost always exists and plays a key role in either stabilizing or destabilizing the interface depending on the relative densities of the uids involved. At a low ow rate the viscous forces within the uids are negligible and the competition between gravity and capillary forces takes control of the displacement. is regime has been studied by (B ., ; A ., ) where gravity : e celebrated "Sa mann-Taylor" ngering observed at the interface between air (top) and glycerine (bo om), showing the development of instabilities and the inhibiting e ect that the fastest growing nger has on its neighbors (S T , ). In the rst image, the velocity of the interface was V = 0.1cm.s -1 and the critical wave length is c r = 1.2cm. e average wave length of the disturbance observed at initiation was 2.2cm which is in the expected range. Note that this experiment was done with the motion of the interface along the direction of gravity.

was found to provide a stabilizing e ect on the drainage front when a denser uid is displaced by a lighter one from the top. On the other hand at a higher ow rate when capillary forces are negligible, the competition between gravity and viscous forces has been studied early on both theoretically and experimentally (in an analogous Hele-Shaw ow context) in the seminal work S T ( ). As an extension to the work of T (

) concerning the condition for instability of a moving interface between two uids of contrasting densities and negligible viscosities in a Hele-Shaw cell, S T ( ) has derived analytically a similar condition applicable to viscous uids in a porous medium: "[. . . ]. When two superposed uids of di erent viscosities are forced through a porous medium in a direction perpendicular to their interface, this surface is stable or unstable to small deviations according as the direction of motion is directed from the more viscous to the less viscous uid or vice versa, whatever the relative densities of the uids, provided that the velocity is su ciently large. "

-S T ( )

For any generic disturbance applied to an interface moving vertically upwards, the above implication was done based on the following analytical condition, stable if :

µ i k i - µ d k d V + (ρ i -ρ d ) > 0, unstable if : µ i k i - µ d k d V + (ρ i -ρ d ) < 0.
( . )

V is the velocity of the interface, the acceleration due to gravity and ρ is the density of the uid under consideration. Note that in the above analysis the invading uid is below the defending uid, in other words the direction of motion of the interface is against that of gravity. Also, no direct consideration has been done as to the we ing properties of the uids involved, to the solid surface and the e ect of surface tension between the uids. However, S T ( ) did indicate, with suggestion from Dr.Chuoke R.L., that the e ect of surface tension, T , is to limit the range of wavelengths of the disturbances to which the interface is unstable. Including surface tension, a linear stability analysis allowed to identify a critical wave length, cr , given by, ) proposed the existence of a macroscopic surface tension, akin to T in Eq.( . ), within porous media ow and thus having the same e ect as in Hele-Shaw ow. Later P F ( ) extended such result to three dimensions. Remark: While Hele-Shaw ow is mathematically akin to two-dimensional ow in porous media, there exist striking di erences especially in the transition region from one uid to another and in the morphology of the ngering instabilities (C W , ; S ., ). In porous media, the transition region from the part of the domain fully saturated with one uid, to the one which is fully saturated with the second uid can be described as a mixing region where both uids partially saturate the pores. On the other hand this transition is sharp in immiscible Hele-Shaw ow. is di erence is primarily due to existence of several internal length scales in the case of porous media (S ., ; H , ). erefore, the viscous ngering instabilities observed in porous media ow have a fractal structure rather than smooth formations as in the case of Hele-Shaw ow. As mentioned earlier, in the derivation of Eq.( . ) gravity was considered acting opposite to the direction of motion of the uid interface. e analysis of the converse case amounts only to changing the sign of acceleration due to gravity giving, stable if :

cr = 2π √ Tb 12 (µ d -µ i ) + b 2 (ρ d -ρ i ) -1 2 , ( . )
µ i k i - µ d k d V + (ρ d -ρ i ) > 0, unstable if : µ i k i - µ d k d V + (ρ d -ρ i ) < 0. ( . )
In fact, more recent works such as M´ . ( ); L . (

), have shown experimentally that gravity, when acting in the direction of motion of the uid interface, helps to stabilize the invasion of a less dense non-we ing uid (ρ i < ρ d ) into a porous medium that is saturated with more viscous (µ i < µ d ) and denser we ing uid. See Fig. ( . ).

Gravity driven in ltration

On the contrary, when a denser (ρ i > ρ d ) and more viscous (µ i > µ d ) we ing uid displaces a resident non-we ing uid, the viscosity contrast is not in favor of destabilization, Eq.( . ), and gravity is expected to assist in destabilizing the ow. is la er regime was the focus and a well-engaged topic of research in soil hydrology since the s. One can look into D C ( ); X ( ) among others for a historical review of gravity driven in ltration of water into soil.

Experimental evidence in this context suggested that ngering type instabilities that occur during in ltration have a non-monotonic pro le along their length, with their tips having higher water content compared to their tails. is was quite evident in the light transmission experiments by G . ( a), Fig. .), due to their absence of macroscopic ngering, allowed to justify the presence of saturation overshoot during in ltration in the range of uxes for which the two-dimensional experiments were found to exhibit ngering. is lead to an e ort of evaluating the new models for their one-dimensional solution pro les in order to predict their ability to describe ngering in higherdimensions D C ( ). On the other hand, this technique could also potentially be useful to determine if ngering could be expected for in ltration into certain types of sands that would otherwise have nger widths limited by the experimental chamber sizes. Only caveat being the edge e ect (D C

., ) that is inevitable due to the contrast between porosities of the bulk and at the edge. Also, it is worth to point out at this moment that, to our knowledge, such quasi one-dimensional experiments

were not yet a empted in the case of drainage. Apart from solution pro les along the direction of in ltration, a feature that is of practical importance is the nger width. is is because instabilities with thinner ngers tends to reduce the e ciency (say of ltration) by completely bypassing large areas of dry soil. Many numerical models intend to exhibit the experimentally observed nger width predictions and their relation with applied in ltration ux. While no comprehensive experimental work exists studying the nger width sensitivity to a large range of ux values, D C ( ) managed to consolidate results from independent two and three dimensional experimental campaigns with appropriate scaling. See Fig. ( . ). Two limits of stable ow were made apparent. At higher uxes, as the ux increased, tending to the saturated conductivity, the nger width increased. In the two-dimensional slab experiments by G . ( b), Fig. ( . ), this increase in the nger width tended to the transverse size chamber. On the other hand, at lower uxes, as the ux reduced the nger width increased as well, tending to stabilized ow. In essence unstable ow is restricted to a certain range of intermediate uxes. is particular feature of the in ltration into dry soil was not predicted by the analytical results according to D C ( ). However, the ability to predict this behavior, e experiments were conducted at increasing ratio, R s , of ux to saturated conductivity ux rates. R s = (0.012, 0.041, 0.088, 0.32, 0.82) to be read from le to right in the gure.

in tandem with the one-dimensional saturation overshoot feature, could be a useful measure to evaluate numerical models. is perspective gains strength due to the correlation between the absence of one-dimensional saturation overshoot at the same extremes of applied ux at which the nger widths tend to be larger, implying stabilization of the front.

. Mode-I fracture openings are of primary concern in geomaterials when it comes to preferential pathways generated due to loss of mechanical integrity of the porous skeleton. Such fractures are prevalent in porous media either saturated or partially saturated. Typical scenarios involving soils are, desiccation cracks due to drying induced evaporation, pressure driven hydraulic fracturing of intact formations, emergence of complex sub-vertical hydrate veins in deep-water marine sediments among others. While the consequences are similar, each scenario has a di erent underlying mechanism that initiates/ propagates the fracture.

In metals owing to their obvious non-granular nature and in cemented materials owing to their cementation, mode-I fracture opening is classically treated as caused by the tensile nature of stresses in the bulk. Typical mode-I bri le fracture treated à la Gri th involves inherent defects on the material surface or the bulk at which fractures initiate and propagate, driven by the release of stored elastic energy which is of a tensile nature.

Unlike metals and cemented materials, soils are inherently granular and the resistance o ered by them is due to their e ective stress dependent cohesion-less frictional behavior. Such 'toughness' of soils has its origins in the con ning stress, inter-granular adhesion due to uid lms in the case of partial saturation, inter-granular friction due to grain texture and packing con guration. Under such conditions, full saturation would mean that there are no uid lms either, that could have otherwise provided some adhesion. So, the sources of resistance against fracturing are the intergranular friction and the con ning stress, σ c . Such is the case of saturated so marine sediments, brine-saturated sub-surface formations and surface soil layers following rainfall/ ooding among others.

When the invading uid is the same as, or miscible with, the one already saturating the porous skeleton, the mechanism behind initiation of fractures at pre-existing defects is the same as that of classical hydraulic fracture. e uid tends to leak-o into the bulk along the normal to the faces of the defect and in proportion to the intrinsic permeability of the porous skeleton. Seepage induced drag forces thus develop within the soil mass around the initial defect. See Fig. ( . ). For fracture to initiate and propagate these drag forces, that are proportional to the di erence in uid pressures within the fracture and the far-eld, must exceed the e ective compressive stresses in the far-eld due to the con nement.

. . Drainage induced fracturing

If the invading uid is non-we ing to the pore walls with respect to the saturating uid, then the mechanism behind a mode-I fracture initiation under compressive e ective stress conditions has been a topic of debate and investigations since the last decade. Indeed, this would be an entry into the realm of partial saturation and a uid-uid displacement that is drainage. A typical example of this scenario is that of desiccation cracks in drying soils.

Desiccation is a process of evaporation induced at the surface of soil that is saturated typically with water. As the surface water evaporates, the water level reaches the soil surface and the air (or water-vapour) resists to invade the surface pores exposed to the exterior. However, air being non-we ing can invade the largest pores only once the capillary pressure di erence between air and water exceeds the entry pressure/ critical capillary pressure. Until such a time the soil mass experiences shrinkage induced by the increasing suction. Once the capillary pressure a ains the entry pressure, air starts to invade the pore spaces. is is when mode-I desiccation cracks seem to be initiated and eventually form networks as observed in various experimental campaigns, notably that of N K ( ); P . One prevailing point of view is that, these fractures are generated under the in uence of tensile forces that are induced either due to boundary conditions that restrain the shrinkage of the whole soil mass (in shallow samples) or due to di erential shrinkage along its depth (in thicker/deeper samples). Relatively recent modeling a empts (C ., ; C S , ; H S , ; H ., ) based on a continuum description of fracture have incorporated such a perspective.

An alternative phenomenon has been championed in the works of S S ( , a,b) and that of J J ( ), which is compatible with the earlier mentioned fundamental behavior of soils, that is to exhibit toughness only under compressive e ective stresses. In this case fracturing is thought to be caused by the forced invasion of the porous network by either uids: immiscible (for instance desiccation) or miscible (earlier explained hydraulic fracturing). In the experiments by S S ( ), a water saturated slurry of clay in a cylindrical chamber has been subjected to pressure by a superposed uid which is either immiscible oil or miscible water. See Fig. ( . ). e miscible case has been setup with an intentional defect that looks like hole. e soil mass is drained in both cases from the bo om.

is is a scenario which unambiguously generates compressive e ective stresses within the con ned sample with the ratio of horizontal to vertical e ective stresses less than 1 that is typical to granular media. As the superposed uid pressure is increased, mode-I fracture initiation has been observed in both cases, starting at surface defects and propagating vertically downward and laterally within a plane normal to the minor principal e ective stress. is de nitely cannot be explained through a tensile stress driven mechanism. In S S ( , a,b) this has been explained through a particle level phenomenology at a supposed surface defect, supported by representative numerical simulations. Numerical simulations for the immiscible case were conducted on a Gaussian defect in a twodimensional domain using a Modi ed Cam Clay model assuming zero soilcohesion and an increasing uniform internal suction, u c (t), applied on the whole domain. e simulations revealed that while consolidation takes place everywhere within the sample accompanied by a reduction in the void ratio, e(x, t), the tip of the defect experiences a lower rate of reduction in void ratio. us as time progresses, in the far-eld the path followed in the space of (u c , e) is that of the normal consolidation line (NCL), whereas the tip deviates from this path thus obtaining a relatively higher void ratio. is higher void ratio corresponds to a lower air-entry suction and thus one comes to a conclusion that air invades rst at the tip of the defect and propagates owing to the same mechanism. e e ective stress, while being lower at the tip, was always found to be in compression throughout the computational domain.

Dependency on grain size and con ning stress A similar explanation of the non-we ing uid invasion has been arrived at by J J ( ) employing discrete element method (DEM) simulations to replicate buoyant methane gas invading marine sediments from below. Here the focus was on the e ect of grain size on the invasion mechanism. Two types of gas invasion into initially saturated granular porous media have been classi ed: classical capillary invasion and invasion by creation of opening mode fracture/s. When the grain size is large, the invasion seem to follow the classical pore-by-pore invasion where the invading gas pressure exceeds the capillary pressure limit at the pore throats which could be su ciently described by invasion percolation theory (L , ). However, when the grain size is reduced, the capillary pressure limit for gas entry at the pore throats is much higher and the increasing invasion gas pressure reaches a fracturing limit whereby the DEM network no longer behaves as a rigid medium. Fractures were found to initiate at the invading boundary and propagate vertically into the network along the direction normal to that of the minor principal e ective stress similar to the experiments in S S ( ). Accordingly a local condition for fracture opening has been identi ed in J J ( ), on the basis of classical fracture mechanics, requiring the gas pressure to exceed the e lateral boundaries were xed and the grain pack is compacted vertically under constant porepressure until vertical e ective stress of 3Mpa is reached. en an incrementally increasing gas pressure is applied at the bottom central pore. Maroon lines depict the compression at the inter-granular contacts, the blue dots are the pores fully occupied by gas and the pink lines are where the intergranular cohesive bonds are broken. (a) Classical capillary invasion for r min = 50µm at capillary pressure of 6KPa. (b) Fracture opening for r min = 0.1µm at capillary pressure of 2.55MPa. minimum compressive stress at the far-eld and the inter-granular adhesion due to uid lms.

H . (

) has carried out an experimental investigation by injecting air into a thin bed of water saturated glass bead pack con ned within a shallow cylindrical chamber with variable weight on top. E ect of the weight that generates a con ning stress, the diameter of the glass beads and the injection ux have been studied by varying them. Apart from the classical ngering transition from capillary to viscous regime with an increasing injection ux, a transition from ngering to fracturing has been observed with reducing con ning stress and as well with reducing bead diameter. See Fig. ( . ). Notably, the la er transition did not seem to be a ected by the varying injection ux.

.

In the case of bi-phasic ow in rigid porous media, even though there is a clear fractal nature to the invasion morphology, at a macroscopic continuum scale destabilization of the invasion seem to be analogous to that of Hele-Shaw ow.

is is further corroborated by the ability of the simplistic stability analysis by S T ( ) being able to qualitatively explain the various regimes of destabilization in the presence of gravity. Moreover, the quasi one-dimensional gravity driven in ltration pro les potentially provide a basis for evaluating any modeling techniques, even before demonstrating their capability to describe ngering in higher dimensions. e prevailing explanation for appearance of opening mode fractures in soils during drainage has been the build up of tensile stresses. e available evidence, both experimental and pore-scale simulations, unequivocally supports an alternate phenomenon, that is of local grain rearrangement caused by the invading capillary lm. Unlike the prevailing explanation, this alternate rationale is inline with the fundamental behavior of granular media like soils, that is an e ective stress dependent frictional strength.

is and the state of current advances within continuum description of fractures, provide suitable grounds for a thorough investigation at a continuum scale by isolating capillarity as the driving mechanism for fracture initiation/propagation. is chapter serves as a recall for two mathematical descriptions of partially saturated poromechanics.

e rst one is the classical approach in which starting from the principles of conservation and the classical constitutive laws, the governing equations are derived. A partial linearization of the governing equations is arrived at, under the hypothesis of small perturbations applicable to poroelasticity. In the special case of a non-deformable porous network, Richards' equation describing the evolution of uids is derived. A phase eld approach is then introduced assuming the presence of two uids as equivalent to a single non-uniform uid fully saturating the pore spaces. While the classical formulation of unsaturated poromechanics is only brie y summarized listing the main hypotheses of the model and the governing equations of the problem, a more detailed deduction of the phase eld approach is discussed. Speci cally the thermodynamic restrictions imposed by an extended form of the Clausius-Duhem inequality to derive admissible constitutive laws for generalized stresses and average uid ux are detailed. Under such a framework the governing equations are derived and an extension to Richards' equation in the non-deformable case is arrived at.

. e following deductions in this section are meant to be a recall to partial saturation as treated in C ( ) using a classical approach. Accordingly a macroscopic treatment is adopted considering the porous medium as a homogenized continuum whose representative elementary volume (REV) is assumed to account for the phase-averaged quantities such as the porosity, saturation degrees, densities, etc.

We consider the pore space to be saturated by two uid phases, one we ing (subscript 'w') and another non-we ing (subscript 'nw') to a deformable skeleton representing the solid phase (subscript 's'). e particular simpli cations that could be obtained by choosing the we ing uid to be water and the non-we ing uid to be air that is typical to soil hydrology are detailed at a later moment in Sec.( . . ). Moreover the following terminology holds for the entirety of this work. A porous medium is the ensemble of the solid skeleton and the uids. A porous solid is a so-called we ed solid skeleton, with a thin layer of immobile we ing uid a ached to its pore walls. Solid matrix is the fabric of the skeleton, which in the current case constitutes the grains of soil.

e primary characteristic of a porous skeleton is its porosity. e Eulerian porosity, n, is de ned as the volume fraction of the pore space within an elementary unit volume of the porous medium, dΩ t at the current time t. Accordingly, partial porosities, n α , α ∈ {w, nw }, distinguish the volume fraction of the pore space occupied by each uid within dΩ t such that n = n w + n nw . e Lagrangian counter parts ϕ α , ϕ are de ned such that n α dΩ t = ϕ α dΩ 0 s and ϕ = ϕ w + ϕ nw with dΩ 0 s representing the unit volume of the porous medium in the reference con guration of the solid skeleton.

Similarly, another useful quantity to de ne is the volume fraction of each uid within the available pore space of dΩ t . is quantity is the saturation degree denoted S α such that 1 = S w + S nw .

If ρ c , c ∈ {s, w, nw } denote the intrinsic mass densities of each phase, then their apparent Eulerian mass densities are given by ρ α n α = ρ α nS α for the uids and ρ s (1 -n) for the solid. e Lagrangian counterparts to the apparent mass densities are given by m α = ρ α ϕ α = ρ α ϕS α for the uids and m s = ρ s (1 -ϕ) for the solid.

For most Eulerian quantities appropriate Lagrangian counterparts exist owing to the deformation of the porous skeleton. Eventually we introduce the hypothesis of small perturbations (HSP) that lets us work under in nitesimal deformations thus unifying both frameworks.

. . Kinematics & preliminaries

Let Ω 0 s and Ω t indicate respectively the reference and current con gurations occupied by the porous medium at time t = 0 and t > 0. Since the porous medium is considered as a homogenized continuum, each spatial position x ∈ Ω t at time t is simultaneously occupied by the solid and the uid particles whose individual reference positions are given by X c , c ∈ {s, w, nw }. e individual deformation maps are given by x = x c (X c , t) and their inverse mappings, with an abuse of notation, by X c = X c (x, t). It is to be noted that it is customary in poromechanics to transport the governing equations of the uids and the solid from their common current con guration, Ω t , to the reference con guration of the solid skeleton, Ω 0 s , so as to allow the usual continuum mechanics approach to constitutive characterization of the solid skeleton.

e time derivatives giving velocity elds of the solid and uid particles identi ed in their current con gurations are,

c (x, t) = ∂x c (X c , t) ∂t X c =X c (x,t )
, c ∈ {s, w, nw }.

( . )

Consequently the relative velocity of the uid particles, X α , w.r.t a solid particle, X s , at every current position x = x s (X s , t) can be obtained as,

w α (x, t) = α (x, t) -s (x, t), α ∈ {w, nw }. ( . )
e Lagrangian gradient w.r.t the solid reference con guration is denoted ∇ s whereas the Eulerian gradient in the current con guration is denoted ∇. Analogously, the divergence operators are denoted ∇ s • and ∇•. With the above notations in place, the deformation gradients of the earlier de ned maps are given by,

F s (X s , t) = ∇ s x s (X s , t), F α (X α , t) = ∇ s x α (X α , t), α ∈ {w, nw }.
( . )

e following transport formulas hold for elementary oriented surfaces, ) and for elementary volumes,

nda = s F -T s Ns dA s , ( .
dΩ t = s dΩ 0 s , ( . )
where s = det F s . nda and dΩ t are the current elementary oriented surface and elementary volume corresponding to the reference counterparts Ns dA s and dΩ 0 s respectively. As usual in continuum mechanics, the displacement of the solid skeleton is measured as u(X s , t) = x s (X s , t) -X s .

e changes in the lengths of material vectors of the solid skeleton and their included angles are measured by the Green-Lagrange strain tensor given by,

E = 1 2 F T s • F s -I = 1 2 (C -I ) = 1 2 ∇ s u + ∇ T s u + ∇ T s u • ∇ s u , ( .
)
where I represents the second-order identity tensor and C is the Cauchy-Green strain tensor.

Particle derivatives

Particle derivative of any di erentiable Eulerian eld G w.r.t the solid or the uid particles is given by,

d c G dt = ∂G ∂t + ∇G • c , c ∈ {s, w, nw }. ( . )
Similarly, the particle derivative applied to a volume integral of any di erentiable eld G is given by,

d c dt ∫ Ω t GdΩ t = ∫ Ω t d c G dt + G∇ • c dΩ t , = ∫ Ω t ∂G ∂t + ∇ • (G c ) dΩ t , = ∫ Ω t ∂G ∂t dΩ t + ∫ ∂Ω t G c • ndS t , c ∈ {s, w, nw }, ( . )
where for the last step the divergence theorem is employed. ∂Ω t is the boundary of Ω t where the outward unit normal is n.

For any Lagrangian eld, G(X s , t), the particle derivative w.r.t the solid, d s G/dt, is further denoted by d G/dt for convenience.

. . Balance of Mass

Following C ( ), in the absence of any sources/ sinks within the bulk and any surface supplies, the Eulerian mass balances applicable everywhere in Ω t for each uid phase α ∈ {w, nw }, and the solid phase in their di erential forms read,

∂(ρ α nS α ) ∂t + ∇ • (ρ α nS α α ) = 0, ∂(ρ s (1 -n)) ∂t + ∇ • (ρ s (1 -n) s ) = 0.
( . )

As a general practice in poromechanics the time derivative in uid mass balance is rewri en w.r.t the solid reference con guration using Eq.( . ) as follows,

d s (ρ α nS α ) dt + ρ α nS α (∇ • s ) + ∇ • (ρ α nS α w α ) = 0. ( . )
With the de nition of the Eulerian relative uid mass ow vector w α (x, t) = ρ α nS α w α = ρ α nS α ( αs ) and the Lagrangian relative uid mass ow vector, M α , such that w α • nda = M α • N s dA s , one can obtain using the transport formulae, the pull-back of the balance of uid masses into the solid reference con guration,

dm α dt + ∇ s • M α = 0, with M α = m α F -1 s • ( α -s ). ( . )
On the other hand, solid mass balance can be rewri en from the integral form of Eq.( . (b)) and then integrating it,

ρ s (1 -n)dΩ t = ρ 0 s (1 -ϕ 0 )dΩ 0 s , ( . )
where ρ 0 s and ϕ 0 are the initial mass density and initial porosity of the solid skeleton. It is apparent that the current density of the solid skeleton is related to the current porosity through an algebraic relation knowing s , ) . .

ρ s = ρ 0 s (1 -ϕ 0 ) s (1 -n) . ( .

Balance of Momentum

In a Galilean reference frame, the instantaneous momentum balance with respect to all the ma er included in the porous domain, Ω t dictates that the sum of variations in time of momenta of each constituent, c ∈ {s, w, nw }, should be equal to the sum of surface and internal supplies. is reads in an integral form as follows for linear momentum, ) with α ∈ {w, nw }. f t and f b are respectively the surface traction acting on the boundary and the body force acting within the domain of the porous medium as a whole. ρ = ρ s (1 -n) + ρ w n w + ρ nw n nw , is the overall mass density of the porous medium. Under a quasi-static assumption, neglecting the dynamic terms in the equation above leads to,

d s dt ∫ Ω t ρ s (1 -n) s dΩ t + d α dt ∫ Ω t ρ α nS α α dΩ t = ∫ Ω t ρ f b dΩ t + ∫ ∂Ω t f t dS t , ( .
0 = ∫ Ω t ρ f b dΩ t + ∫ ∂Ω t f t dS t . ( . )
Classical Cauchy stress theorem (C , ) linearly relates the surface traction to the outward unit normal at the boundary where it is applied through the identi cation of a Cauchy stress tensor, σ (x, t), for the overall porous medium, such that f t = σ • n. Using this in Eq.( . ) results in,

0 = ∫ Ω t (ρ f b + ∇ • σ ) dΩ t .
( . )

Whereas, the balance of angular momentum as usual in continuum mechanics would amount to the symmetry of the tensor σ (x, t).

. . Hypothesis of small perturbations (HSP)

e hypothesis of small perturbations in the context of poromechanics amounts to a set of assumptions that allow a partial linearization of the highly non-linear problem. ese are listed below in no particular order:

• Hypothesis of in nitesimal transformations:

∇u 1. ( . )
According to which the Lagrangian and Eulerian approaches coincide up to a rst-order approximation. It also implies in nitesimal deformations giving E 1. As a consequence the Green-Lagrange strain tensor, E, can be approximated by the linearized strain tensor, ε, given by,

E(u) ≈ ε(u) = 1 2 ∇u + ∇ T u . ( . )
• Hypothesis of small displacements: is states that the displacements of the solid skeleton are much smaller than the characteristic length scale, L, of the structure as a whole,

u/L 1. ( . )
is amounts to a merger of the reference and the current con gurations Ω 0 s ≡ Ω t = Ω and their respective boundaries ∂Ω 0 s ≡ ∂Ω t = ∂Ω. Such a merger is only as far as the spatial arguments of unknown elds are concerned, i.e., x X whenever they appear as arguments of unknown functions. As result of the above two hypotheses, the operator ∇ can be used directly without speci c reference to the con guration, ∇ ≡ ∇ s .

• Hypothesis of small variations of Lagrangian porosity: is amounts to,

|(ϕ -ϕ 0 )/ϕ 0 | 1. ( . )
It is to be noted that this is not equivalent to a non-deformable porous skeleton, ϕ = ϕ 0 , which is treated at a later point.

• Hypothesis of small variations of uid mass densities: is amounts to,

|(ρ α -ρ 0 α )/ρ 0 α | 1, α ∈ {w, nw }. ( . )
Superseding this hypothesis, uids in liquid state are assumed incompressible and uniform in the current work, such that ρ α (x, t) = ρ 0 α .

. . Classical constitutive relations

With a choice of minimum primary unknowns, the balance laws are supplemented by constitutive relations that are derived employing the hypothesis of energy separation (see Sec.( . . ) in C ( )). According to this, the individual contributions of the free energy density of the porous skeleton, ψ s (ε, ϕ), and the interfacial energy density, U (S w ), are separated from within the free energy density of the porous solid, Ψ s (ε, ϕ, S w ),

Ψ s (ε, ϕ, S w ) = ψ s (ε, ϕ) + ϕU (S w ).
( . )

In doing so, U (S w ) is understood as encompassing both the solid-uid and uid-uid interfacial contributions. Starting from ψ s (ε, ϕ), a er a Legendre transformation involving the conjugate pair (ϕ, π ), an alternative state function describing a partially saturated deformable porous skeleton is obtained as G s (ε, π ) where π is an equivalent pore pressure described below.

e corresponding state equations are the constitutive laws that account for the tensile e ect of the interfaces on the pore walls. e material properties which specify the linear poroelastic model are typically identi ed through experimental investigations.

Without resorting to exhaustive derivations, below we recall brie y these relations under an overarching assumption of isotropic linear elastic homogeneous porous medium undergoing isothermal transformations. e primary unknowns chosen for this purpose are the displacement of the solid skeleton, u, and the uid pressures, p α , α ∈ {w, nw }.

Fluid density e state equations for the uid phases, α ∈ {w, nw }, relate their intrinsic densities to their averaged uid pressures p α assuming simplistic barotropic uid phases,

ρ α = ρ α (p α ). ( . )
Under the HSP in-compressible uids are assumed and the above simply reduces to ρ α = ρ 0 α , where ρ 0 α is the corresponding initial intrinsic uid mass density.

Capillary pressure & interfacial energy

At a pore-scale, the menisci formed between the we ing and non-we ing uids have certain curvature depending on the surface tension between the two uids, their we ability w.r.t the pore walls and the geometry of the pores. e pressure di erence between the two uids, that is the porescale capillary pressure, governs the volume fractions of those uids within the pore space. Classical poromechanical practice involves assuming a simplistic geometry of the pore space and specifying that the work done by the pore-scale capillary pressure in causing in nitesimal variation in the volume content of the we ing uid occupying the pore, is equivalent to the in nitesimal variation of interfacial energy (C , ). us a continuum scale constitutive relation is determined between the interfacial energy, U (S w ), and the saturation degree of the we ing uid, S w . See Fig. ( . ).

is relation is called the retention relation or the retention curve, that at equilibrium allows for the retention of a given volume of we ing uid at a prescribed pressure di erence between the two phases, (p nw -p w ) = p c (S w ), within a porous skeleton. In this sense, it accounts not only for the surface tension between the two uids but also for the retention e ect provided by the pore walls due their texture.

e functional form of p c (S w ) usually is empirical ( G , ; B C , ) in order to t the experimental data and can generally be hysteretic. In the current study the widely used van Genuchten form ( G , ),

p c (S w ) = - ∂U ∂S w = π 0 S w -S r es w 1 -S r es w -1 m -1 1-m , ( .
)
is adopted where, π 0 scales as ρ w c . c [L] is the porous media counterpart of capillary length scale that relates intensity of capillary forces in the porous medium to gravity. S r es w is the saturation degree of the residual we ing uid that is immobile and is always a ached to the pore walls. Also, as can be observed from Fig. ( . ), p c (S w ) is a non-negative function of S w ∈ [S r es w , 1] that monotonically decreases from +∞ at S w = S r es w to 0 at S w = 1. So, U (S w ) has a unique minimum within S w ∈ [S r es w , 1] at S w = 1. It is to be noted that since uid pressures, p α , are chosen as the primary unknowns, saturation degree of the we ing uid can be obtained by the inverse relation associated to Eq.( . ), ) wri en with an abuse of notation. e saturation degree of the non-we ing uid can then be obtained as, S nw = 1 -S w . Also, as a cautionary note, the above relation assumes a homogeneous retention property, whereas spatial variations of porosities due to skeleton deformation could invalidate such an assumption. e validity of uniform retention properties is related to the HSP and is applicable only to the framework of linear poroelasticity.

S w (p w , p nw ) = p -1 c (p nw -p w ), ( .

Equivalent pore pressure & average fluid pressure

In full saturation where a single uid completely saturates the pore spaces, the pore pressure is equivalent to the uid pressure. In partial saturation however, due to the presence of two uids phases superposed with the solid phase in the current con guration, a need for an equivalent pore pressure, π , arises. is equivalent pore pressure should properly take into account the contributions of uid pressures of each uid phase supposedly acting in unison on the pore walls.

While various possibilities for such an equivalent pore pressure were proposed over the years, we follow here the thermodynamically justi ed version by C ( ),

π (p w , p nw ) = p w S w + p nw S nw -U (S w ). ( . )
Such a de nition accounts not only the e ect of an average uid pressure , p * = p w S w + p nw S nw , but also for the tensile e ect on the pore walls due to the capillary interface.

E ective stress

Since the early works of Terzaghi (T , ) in soil mechanics and of Biot (B , ) in rock mechanics the concept of an e ective stress is introduced that is assumed responsible for the entirety of the deformations of the solid skeleton, while taking into account the contribution due to the uid pressure in the case of full saturation. An extension to the case of partial saturation for in-compressible solid matrix proposed by Bishop (B , ) takes into account the concept of equivalent pore pressure, π . While originally proposed phenomenologically, the e ective stress, σ e , can be derived also through thermodynamic considerations as done in C (

). Following the notation in the la er it is de ned through,

σ -σ 0 = σ e -b(π -π 0 )I, ( . )
where σ 0 is the total pre-stress and π 0 is the equivalent pre-pressure if any. b is the so-called Biot's coe cient which in the case of homogeneous isotropic solid skeleton can be obtained (B W , ) from the bulk modulus of the empty skeleton K and the bulk modulus of the solid matrix

K s , b = 1 - K K s . ( . )
Assuming an isotropic linear poroelastic skeleton results in the e ective stress being related to the linearized strain tensor, ε(u), through the fourth order sti ness tensor C as follows, . ) which is equivalent to an isotropic Hooke's law for classical materials. ϵ 0 is the pre-strain if any. Due to the said assumption the sti ness tensor can be fully determined with the knowledge at least two material constants related to the empty porous skeleton, say the rst, λ = K -2 3 µ, and the second µ Lamé parameters, µ being the corresponding shear modulus.

σ e = C . . (ε -ε 0 ), (
us the following classical relation holds,

σ e = λϵ(u)I + 2µε(u), ( . )
with ϵ(u) = tr(ε(u)).

Porosity

Again adopting an isotropic linear poroelastic solid, change in the total Lagrangian porosity, ∆ϕ, is obtained as a sum of contributions due to variations in solid strain and in equivalent pore pressure. is is given by the relation,

∆ϕ = ϕ(u, p w , p nw ) -ϕ 0 = b(ϵ -ϵ 0 ) + 1 N (π -π 0 ), ( .
)
where N is the so-called Biot's modulus that can be related to the properties of the macroscopic skeleton and of the matrix through,

1 N = b -ϕ 0 K s . ( . )
Remark: In soils the individual grains are typically in-compressible implying K s goes to in nity, b ≈ 1 and 1/N ≈ 0. Consequently, changes in porosity are entirely owing to skeleton strains, ∆ϕ ≈ ϵ(u). e e ective stress, σ e , reduces to the Bishop's e ective stress, σ + πI . In the case of full saturation this is equivalent to the simpli cation of Biot's e ective stress to the celebrated Terzaghi's e ective stress.

Darcy's law e biphasic uid ow accounting for dissipations of the uid phases is described by an extension of Darcy ow to the partially saturated conditions.

is is a relation between the volumetric ux of each uid phase w.r.t the solid skeleton and the driving force which usually is the gradient of the uid pressure of the respective phases and gravity (if considered). is is given for the each uid phase as, ) where κ and η α are respectively the intrinsic permeability of the porous skeleton and the dynamic viscosity of the uid phase α. is the acceleration due to gravity. V α is the so-called ltration vector that represents the volumetric ux of each phase relative to the solid skeleton. k α (S α ) represents the relative permeability function for each phase which is a non-linear and typically empirical relation that modulates the ow of the particular uid phase in accordance to its local saturation degree. In the current work we adopt empirical relations given by the van Genuchten model ( G , ; L ., ),

w α ρ α = V α = nS α ( α -s ) = κ η α k α (S α ) (-∇p α + ρ α ) , ( .
k w (S w ) = S w 1 -1 -S w 1 m m 2 , ( . ) k nw (S nw ) = S nw 1 -(1 -S nw ) 1 m 2m . ( . )
e above relations are smooth functions, see Fig. ( . ), that behave such that the more saturated the porous network is with a particular phase α the closer to one the relative permeability, k α (S α ), of that phase will be. It is to be noted that the function k w (S w ) is typically purely convex within S w ∈ [0, 1].

. . Summary of Governing equations under HSP

Under the simpli cations provided by HSP, the governing equations are linearized to an extent. e residual non-linearity comes from the constitutive relations related to the uids and the unsteady term within the balance of mass.

• e mass balance of the uid phases, α ∈ {w, nw }, given by Eq.( . ) reduces to,

∂(ϕS α ) ∂t + ∇ • V α = 0, ( . )
with V α given by the Darcy's law with ρ α = ρ 0 α , recapitulated below for convenience,

V α = κ η α k α (S α ) -∇p α + ρ 0 α . ( . )
• e evolution of intrinsic mass density of the solid skeleton given by Eq.( . ) reduces to,

ρ s = ρ 0 s (1 -ϕ 0 ) (1 -ϕ) . ( .
)
Since this is an algebraic relation, ρ s could be derived once ϕ is known from the solution of the problem through Eq.( . ).

• e linear momentum balance or the equilibrium equation given by Eq.( . ) with the help of classical localization arguments reduces to the di erential form,

∇ • σ + ρ f b = 0, ( . )
with appropriate boundary conditions listed further. e total stress, σ , is given by Eq.( . ) and the e ective stress, σ e , by Eq.( . ) which have already been stated under HSP. e overall mass density of the porous medium is now given by ρ = ρ 0 s (1 -ϕ 0 ) + (ρ 0 w S w + ρ 0 nw S nw )ϕ. f b in principle accounts for the ensemble of body forces acting the porous medium which could also account for gravity.

. . Boundary & initial conditions

Appropriate boundary (BC) and initial (IC) conditions are needed to close the system of equations that form a boundary-initial value problem (BIVP).

Typically the ICs relating to the state of the porous medium are assumed to be homogeneous ones. For instance, in soil hydrology one would usually come across initially dry or fully saturated conditions. e initial displacement corresponds to pre-strain of the solid skeleton if any. ese set of ICs read as follows at t = 0 and for any x ∈ Ω, u(x, 0) = u 0 , p α (x, 0) = p 0 α , α ∈ {w, nw }.

( . )

e BCs typically appear as essential (Dirichlet), natural (Neumann) and mixed (Robin) conditions. For the purpose of the current work we only consider the rst two types of BCs. Parts of the full boundary, ∂Ω, where these two types of BCs are applied are designated with a subscript 'D' and 'N ' respectively. e mechanical BCs complimenting the equilibrium, Eq.( . ), of the porous medium are designated with a superscript 's' and those complimenting the balances of uid mass, Eq.( . ), are designated with a superscript 'α' corresponding to the uid phase concerned, α ∈ {w, nw }. Further, we assume that all the BCs are independent of time and are uniformly applied to their corresponding parts of the boundary. ese set of BCs read as follows for all t > 0,

u(x, t) = u D , x ∈ ∂Ω s D , ( . ) σ • n = f N , x ∈ ∂Ω s N ( . ) p α (x, t) = p D α , x ∈ ∂Ω α D , α ∈ {w, nw }, ( . ) V α • n = q N α , x ∈ ∂Ω α N , α ∈ {w, nw }. ( . )
With an abuse of notation n denotes the outward unit normal to the boundary wherever invoked. σ and V α are given by their constitutive relations, Eq.( . ) and Eq.( . ) respectively. u D , f N , p D α and q N α are respectively prescribed values, on the corresponding boundaries, of displacement, surface traction, uid pressures and volumetric uid uxes in the direction of n.

. . Soil hydrology & Richards' equation

Typically in soil hydrology the two uids involved are water and wet-air (or water vapor). Water in its liquid state is considered in-compressible, ρ w (x, t) = ρ 0 w , which is already taken into account until this moment. Air on the other hand is naturally incompatible with such an assumption and rather can be modeled as an ideal gas,

ρ nw (x, t) = p nw (x, t) R T , x ∈ Ω ( . )
where R and T are respectively the speci c gas constant and absolute temperature assumed uniform and constant. As a rst approximation considering air to be highly rare ed implies that its density is vanishing, ρ nw (x, t) ≈ 0 and an assumption of continuous connection of the air phase to the exterior ensures that any gradients in pressure are instantaneously vanishing. In fact, from the above equation one gets p nw (x, t) ≈ 0. Introducing these into the Darcy's law for air phase in Eq.( . ) gives V α ≈ 0. us the unknowns related to the air phase are eliminated and its saturation degree can be obtained trivially, S nw = 1 -S w , once the problem related to the water phase is resolved. Concurrently, the capillary pressure reduces to p c = -p w and the equivalent pore pressure to,

π = p w S w -U (S w ). ( . )
e overall mass density of the porous medium is now approximated by ρ = ρ 0 s (1 -ϕ 0 ) + ρ 0 w S w ϕ. Now coming to the problem related to the water phase, inserting the Darcy's law Eq.( . ) into the balance of mass, Eq.( . ) gives Richards' equation,

∂(ϕS w ) ∂p w ∂p w ∂t + κ η w ∇ • kw (p w ) -∇p w + ρ 0 w = 0, ( . )
typically a ributed to R ( ) but rst published by R (

). In the above writing the dependency on the primary uid unknown, i.e. p w , is explicitly shown and hence this is referred to as the "head" form of Richards' equation. Note that S w (p w ) is obtained through the inverse retention relation, Eq.( . ) and the relative permeability function is rewri en with p w as the argument, hence is denoted kw (p w ).

e Richards' equation can also be wri en in the "mixed water content" form by considering S w as the primary unknown and water pressure as its function, p w (S w ) = -p c (S w ) wherever it appears, leading to,

∂(ϕS w ) ∂t + κ η w ∇ • k w (S w ) ∇p c (S w ) + ρ 0 w = 0. ( . )
In practical applications where soil is made up of layers with di erent properties, the water content among these layers undergoes sharp changes while the pressure remains continuous. is renders the head form more suitable. Whereas, in uniform soils the mixed water content form is more amenable, especially because use of the head form requires high accuracy time integration to avoid mass balance errors (T ., ) due to the high non-linearity involved in the coe cient of the unsteady term, ∂(ϕS α )/∂p w . Either of these forms are solved in tandem with the equilibrium equation Eq.( . ), appropriate ICs and BCs (except those of air phase) to address water movement in deformable soils. In the case of non-deformable skeleton the porosity remains unchanged (ϕ = ϕ 0 = const) and the only problem to be resolved is of the ow of water and Richards' equation reduces to,

ϕ ∂S w ∂t + κ η w ∇ • k w (S w ) ∇p c (S w ) + ρ 0 w = 0. ( . )
. .

Criticism and extensions of Richards' equation

While this approach has been widely accepted, over the years Richards' equation has faced criticisms on several fronts even though it reduces the complexity involved. ese criticisms surround: the rough up-scaling law that is involved to model the retention e ect of the pore walls on the uids (M , ; G H , a), inability to reproduce physically observed ngering phenomenon (E ., ; N .,

) and complexities involved in numerical implementations M . (

); F O ( ). Concerning the la er, one can notice that the two empirical constitutive relations for the relative permeability, k w (S w ), and the capillary pressure, p c (S w ), are highly non-linear and are degenerate. k w (S w ) can behave as a constant close to full saturation or vanish in close to dry conditions. p c (S w ) as well takes arbitrarily small values close to full saturation and behaves asymptotically when tending to residual saturation. In reality these relations can as well be found to be hysteretic. In the case of in ltration in to dry soils two main features have been observed of the we ing phenomenon: sharp we ing front and ngering type instability as the front propagates through the soil mass (see Sec.( . . )). e extremely large spatial gradients of the solution render a high degree of non-linearity in Richards' equation making the design of numerical schemes complicated (M ., ). In terms of stability, the works of E . ( ); N . ( ) revealed that Richards' equation is unconditionally stable against traversal perturbations thus being unable to produce experimentally observed ngering type instabilities.

In addition to the above criticisms, one needs to practice caution when using Richards' equation or its associated assumptions.

e continuous connectivity of air phase to the exterior and the absence of impermeable layers is critical for neglecting the ow of air. Moreover, in the degenerate case of vanishing saturation degree of water, S w → 0, the equivalent pore pressure does not vanish. On the contrary it tends to negative in nity in accordance with the behavior of capillary energy, π = -p c (S w )S w -U (S w ) → -∞, and with respect to the in nitely thin residual pore water lms a ached to the grains. However, the pore pressure within the dry soil should be equal to the atmospheric pressure due to the continuous connectivity of the air phase which is not compatible with the behavior of equivalent pore pressure.

While various extensions have been proposed intending to remedy these inadequacies, we mention here a few of them which focus on enriching the retention relation and enabling the description of ngering phenomenon. H G ( ); G H ( a,b) introduced the speci c interfacial area, which is a cumulative measure of uid interfaces within an REV, as an internal state variable thus providing a corrective term to otherwise coarsely de ned retention relation, p c (S w ). As part of their 'hold-back-pile-up' postulate E G ( ) drew parallels to this corrective term, introducing the so-called 'hypo-di usive' term and showed (E G , ) that if this correction results in a non-monotonic retention relation then non-monotonic overshoot saturation pro les are possible in gravity driven in ltration. E G ( ) also introduced a mixed 'relaxation' term, which is second-order in space and rst-order in time, as an extension to standard Richards' equation ) and showed similar overshoot solutions. C F J ( a) introduced a phase eld model with a fourth-order in space term to the standard Richards' equation producing non-monotonic solution pro les and proposed a relevant heuristic scaling of the fourth-order coe cient that resulted in favorable comparison with some experimental observations. e caveat of this model is the use of a nonphysical 'compressibility' term of an exponential form thus modifying the behavior of the capillary energy. Stability analysis (C F J , b) of this model as well has revealed instabilities when saturation pro les are non-monotonic.

. More recently in the work of S ( ), a framework has been built that is thermodynamically consistent and naturally allows for modeling possible localized deformations of the porous skeleton (V ., ; D V , ; D A ` , ) and its coupling to the uid ow. Concerning the uid description, while the standard retention relation describes the con ning e ect due to the pore walls, an enhanced description of surface tension between the uids is obtained by considering the two immiscible uids as a single non-uniform biphasic uid in the spirit of C H ( , ). In this sense it can be viewed as a phase eld approach to model multi-phase uid ow while accounting for the fact that the porous network, through which the non-uniform uid is owing, provides additional interfacial energy due to the solid-uid interfaces. Such a description of uids is adopted in Chp.( ). Naturally within such a framework, coexistence of isopotential phases is possible thanks to the Van der Waals-like double-well form of the uid energy density. is gives rise to a non-monotonic e ective retention curve when the con ning e ect given by the standard retention curve is added to it.

is is elaborated further in this chapter. Parallels may be drawn to the 'hypo-di usive' term (E G , ; D C ., ) which resulted in a non-monotonic retention relation, but the mixed 'relaxation' term is not employed in the current model. Instead, as usual in Cahn-Hilliard type phase eld models (L T , ; J , ; B L , ; K , ), coarsening of the uid domains and pa ern formation is made possible due to a non-local gradient energy contribution giving rise to an extension to Richards equation which is now fourth-order in space. In this sense parallels may be as well drawn to the higher order model introduced by C F J ( a), which, however, does not account for coexistence of isopotential phases and phase segregation due to its lack of a bulk energy contribution that has a double-well structure.

In the following sections a derivation of this model is done following S ( ), however neglecting the description of a second-gradient solid through a simplifying assumption on the relation between solid and uid double forces.

. . Kinematics & preliminaries

A macroscopic approach is adopted for treating the kinematics of a porous medium, which is understood, similar to classical poromechanics, as a superposition of the solid skeleton and the mixture of constituent uids (water and air). As detailed further, the binary uid mixture is assumed to be a single non-un orm uid in the same sense as Cahn-Hilliard uids, wherein possible spatial variations of an intensive scalar property is expected to describe the distribution of the constituent uids. With such an assumption, the current kinematic description is closer to the classical one describing full saturation of the pore spaces with a single uid (see C ( )). As a consequence and in view of being concise, the kinematics of classical partial saturation as described in Sec.( . . ) is adapted for the following. Speci cally, the two constituents, c ∈ {s, f }, of the porous medium are considered the solid phase (subscript 's') representing the skeleton and the uid phase (subscript 'f ') representing the non-uniform uid. With this, all the developments of Sec.( . . ) can be applied with subscript 'α' replaced with subscript 'f '.

As mentioned earlier, an intensive property is needed to delineate the constituent phases (water and air) within the non-uniform uid,thus giving their spatial distribution. Following S ( ), this property is chosen as the ratio of volume density of the liquid phase (water) to the available porous volume, which is the nothing but the saturation degree of the liquid phase. is is denoted as S r whose notational counterpart in the classical approach is S w . In the spirit of Cahn-Hilliard modeling, S r acts as a phase eld/order parameter distinguishing the two phases. is particular choice puts a physically motivated restriction on the range of values this phase eld parameter can take, as such S r ∈ [0, 1]. S r = 0 represents pure phase of air(non-we ing uid) and S r = 1 represents water(we ing uid). Moreover, it allows to invoke the assumptions (H S , ) that lead to Richards' equation. For an S r ∈ [0, 1], a simple linear interpolation of individual densities of pure phases can be assumed for the apparent mass density, ρ f , of the non-uniform uid (K ,

). A non-uniform uid composed of in-compressible water (ρ w = const) and highly rare ed air, the la er being continuously connected to the atmosphere, implies that pure air phase is of in nite mobility. us any pressure di erences in this phase (S r = 0) are instantaneously vanishing and the pressure itself equates to that of reference atmospheric pressure, assumed 0 for convenience. Concurrently invoking an ideal gas law for air implies that its density, ρ nw , is negligible compared to that of water (ρ w ρ nw ≈ 0) leading to,

ρ f = ρ w S r + ρ nw (1 -S r ) ≈ ρ w S r . ( . )
. .

Balance of Mass

Similar to the developments in Sec.( . . ), the balance of masses wri en for the non-uniform uid and the solid phases would result in,

dm f dt + ∇ s • M = 0, with M = m f F -1 s • ( f -s ) = m f F -1 s • w, ( . )
as the pull-back in the solid reference con guration of the uid mass balance and Eq.( . ) as the relation for the current mass density of the solid.

. . Power of external forces & balance of momentum

According to the classical gradient theories, the power of external forces is wri en as a linear functional of the velocity elds of the non-uniform uid and the solid,

P ex t = c ∫ Ω t b c • c dΩ t + ∫ ∂Ω t t c • c dS t + ∫ ∂Ω t τ c • ∂ c ∂n dS t , ( . )
with c = {s, f }. ∂Ω t is the smooth boundary of Ω t such that the outward normal n to the boundary does not su er jumps. Note that even though the deformations of the solid skeleton are eventually expected to be described using the standard rst gradient formulation, we start the deductions by taking into account the contribution of a solid double force, τ s , in the expression above, in addition to the bulk force, b s , and the surface traction, t s . Whereas a contribution due to the uid double force, τ f , is considered in addition to its bulk force, b f , and surface traction, t f , owing to its gradient nature elaborated further. In line with the extended Cauchy theorem (D 'I S , ) the surface tractions and double forces for each constituent can be related to the outward unit normal through tensorial quantities as follows,

(Σ c -∇ • Π c ) • n -∇ s • (Π c • n) = t c , on ∂Ω t , ( . ) (Π c • n) • n = τ c , on ∂Ω t . ( . )
Σ c and Π c are respectively second-order apparent Cauchy stress and thirdorder apparent hyper-stress tensors per unit volume of the porous medium in the current con guration. Introducing these into Eq.( . ) one gets,

P ex t = c ∫ Ω t b c • c dΩ t + ∫ ∂Ω t [(Σ c -∇ • Π c ) • n] • c dS t + ∫ ∂Ω t (Π c • n) . . ∇ c dS t , = c ∫ Ω t b c • c dΩ t + ∫ Ω t [∇ • (Σ c -∇ • Π c )] • c + Σ c . . ∇ c dΩ t + ∫ Ω t Π c . . . ∇∇ c dΩ t .
( . )

According to the classical poromechanical practice we de ne the overall stress acting on the porous medium as a sum of the individual stresses acting on the solid and the uid, Σ = Σ s + Σ f , similarly for the overall hyper-stress, Π = Π s + Π f , and the overall bulk force as the sum of the individual bulk forces, b = b s +b f . A quasi-static overall balance of momentum of the porous medium taking into account the extended Cauchy theorem reads in its local form,

∇ • (Σ -∇ • Π ) + b = 0, in Ω t .
( . )

At this point we make a couple of simplifying assumptions. Firstly, we assume that the structure of the hyper-stress tensor for both constituents is such that Π c = I ⊗ π c , where π c is equivalent to an hyper-stress vector. Such an assumption amounts to a restriction on the double force such that it works only on di erential elongation in the case of the solid phase and on di erential saturation ratio in the case of the uid. is means that the skew-symmetric couples working on the vorticity on the boundary due to τ c are assumed to be negligible. e second assumption concerns the overall hyper-stress itself, which is vanishing by assuming Π f = -Π s . e e ect of this assumption as seen further is the absence of second-gradient solid contribution in the expression of strain working. Consequently the balance of momentum reduces to the classical one,

∇ • Σ + b = 0, in Ω t . ( . )
Employing the principal of virtual power (P ex t = P int ) and introducing Eq.( . ) into Eq.( . ) we can get, following a trivial manipulation, an expression for the power of internal forces or the strain working as,

P int = ∫ Ω t b f • w + ∇ • Σ f -∇ • Π f • w + Π f . . . ∇∇w + Σ . . ∇ s dΩ t .
( . )

In order to perform a pull-back of the power of internal forces onto the solid reference con guration one needs to introduce appropriate Lagrangian quantities. e overall second Piola-Kirchho stress tensor, S, as a counter part to the overall Cauchy stress, Σ, is obtained such that,

s 2 Σ . . (∇ s + ∇ T s ) = S . . E =⇒ Σ = -1 s F s • S • F T s , ( . )
where ( ) represents time derivative w.r.t to the solid particle xed in the solid reference con guration. e pull-back into the solid reference con guration of the uid bulk forces is obtained as b f = s F T s b f . Further the uid stress tensor can be decomposed into its spherical and deviatoric parts as, Σ f = -p f I + Σ d f . Since the uid hyper-stress vector is de ned such that π f = Π f /3, its pull-back, γ f , is obtained in accordance with π f = -1 s F s • γ f . Introducing these quantities, the Lagrangian writing of P int , is obtained as,

P int 0 = ∫ Ω 0 s b f • M m f -∇ s • s p f M m f + ∇ s • s F -1 s • Σ d f T • F s M m f -∇ s • ∇ s • γ f M m f - γ f s ∇ s • s M m f + S . . E dΩ 0 s .
( . )

. . ermodynamics

As mentioned earlier (Sec.( . . )), the classical approach to unsaturated poromechanics, as detailed in C ( , ), involves a rough upscaling of a simplistic pore-scale mechanism: energy stored due to creation of additional in nitesimal interfacial area within the pore is due to the work done by the capillary pressure di erence causing in nitesimal variations in the volume content of the we ing uid with the pore. Criticism concerning this approach and the particular extensions to remedy this lacking are detailed in Sec.( . . ).

Fluid & interfacial energies

In the current thesis we follow a novel approach proposed by S ( ) that involves a separation of the uid-uid and solid-uid interfacial energies.

e former is taken into account in the macroscopic constitutive prescription by describing the uid mixture (air and water) as a single bi-phasic non-uniform uid analogous to a Cahn-Hilliard (C H , ) like (L T , ; J , ; B L , ; K , ) uid modeling. Accordingly the non-uniform uid's internal energy is assumed to be a sum of local and non-local contributions,

E f = nρ f e f 1 ρ f , s f + κ f (f p ), f p = ∇(nρ f ) • ∇(nρ f ) . ( . )
e f is the uid-speci c internal energy. ρ f e f is assumed to be a double-well potential dependent on the speci c density, ρ f and speci c entropy, s f of the non-uniform uid. κ f (f p ) is the non-local energy dependent on the gradient of the non-uniform uid content, so is responsible for regularization of the non-convex local energy and penalization allowing formation of interfaces. With ψ f as the uid-speci c Helmholtz free energy and according to the Legendre transformation ψ f = e f -T s f , thermodynamic pressure, P, and chemical potential, µ, of the uid can be de ned such that,

P = - ∂e f ∂ 1/ρ f s f =const = - ∂ψ f ∂ 1/ρ f T =const , µ = ∂(ρ f e f ) ∂ρ f T =const , ( . )
where T is the absolute temperature, conjugate to s f .

In view of the constitutive prescription of the uid, Eq.( .), the following expressions for uid stress and hyper-stress are derived,

p f = nP -κ f + 2 1 + 1 trI ∂κ f ∂ f p f p , ( . ) Σ d f = -2 ∂κ f ∂ f p ∇(nρ f ) ⊗ ∇(nρ f ) - 1 trI ∇(nρ f ) • ∇(nρ f )I , ( .
)

π f = -2 ∂κ f ∂ f p (nρ f )∇(nρ f ). ( . )
e local double-well energy, similar to Van der Waals' isotherms in phase transition, is intended to allow for coexistence of the immiscible phases corresponding to isopotential minima. In view of the assumptions leading to the de nition Eq.( . ) for mass density of the non-uniform uid in terms of the saturation degree of the we ing uid, a possible form of the double-well energy potential is chosen as a function of S r ,

ρ f e f = Cγ R S 2 r (1 -S r ) 2 . ( . )
is form has a symmetric double-well structure whose isopotential minima e characteristic radius of channels within which the uid menisci form, R, is an intrinsic property of the porous skeleton and according to L ( ) it scales with κ/ϕ 0 . With Ψ f (S r ) as the free-energy counterpart of the above local internal energy,

µ f (S r ) = ∂Ψ f (S r ) ∂S r T =const = ∂(ρ f e f ) ∂S r T =const = 2Cγ R S r 1 -3S r + 2S 2 r , ( . )
de nes a chemical potential that is identical to ρ w µ. e non-local gradient energy contribution allowing the formation of di use interface between the two phases and providing regularization of the non-convex local energy, is typically assumed to be quadratic in the gradient of uid content, ∇(nS r ), as follows:

κ f (f p ) = C k 2 (∇(nS r ) • ∇(nS r )) . ( . )
e magnitude of C k determines the spatial in uence of this non-local energy. e signi cance and choice of scaling of this term is reserved for later analysis. However, in the case of a homogeneous evolution of S r this term vanishes. Introducing Ψ N L as the free-energy counterpart of κ f (f p ),

µ N L (∇(nS r )) = ∇ • ∂Ψ N L ∂∇(nS r ) T =const = ∇ • ∂κ f ∂∇(nS r ) T =const = -∇ • π f nS r , ( . )
de nes the non-local chemical potential that is wri en here in view of its future use in Sec.( . . ).

Concerning the solid-uid interfacial energy, E sf , a priori no constitutive prescription is done. is interfacial energy is thought of as responsible for the con ning e ect on the non-uniform uid, provided by the pore walls due to their texture. Following the thermodynamic restrictions this contribution is to be found equivalent to the classical interfacial energy, U (S r ).

The rst principle of thermodynamics

Assuming the postulate of local state extended to porous continua (see C ( )) and introducing E s = ρ s (1 -n)e s as the internal energy of the solid skeleton per unit volume of the porous medium, the rst principle of thermodynamics amounts to a balance of energy that reads as,

d s dt ∫ Ω t E s + 1 2 E s f dΩ t + d f dt ∫ Ω t E f + 1 2 E sf dΩ t = P int + Q. ( . )
Note that the solid-uid interfacial energy is assumed to be equally decomposed in such way as to associate half of it to the displacement of solid matrix and the other half to the motion of the non-uniform uid. Assuming no volume heat rate supply, Q denotes the rate of heat supply from the exterior through the surface of Ω t . is can be de ned, employing the ux-density lemma, in terms of a rate of surface ux density, q, supplied through ∂Ω t along its outward unit normal as,

Q = - ∫ ∂Ω t q • n dS t = - ∫ ∂Ω 0 s q • Ns dS 0 s = - ∫ Ω 0 s ∇ s • q dΩ 0 s , ( .
)
where q is introduced as the pull-back into the solid reference con guration of q.

An overall internal energy per unit volume of the porous medium can be de ned as E = E s + E f + E s f . Employing the de nitions of particle derivative of an integral, Eq.( . ), and the relative uid mass ow vector, w(x, t), the following can be derived, ) where in the last step a pull-back into the solid reference con guration is performed with the introduction of corresponding overall internal energy per unit volume of the porous medium such that EdΩ 0 s = EdΩ t . Finally rewriting Eq.( . ) in the solid reference con guration using, Eqs.( . , . & . ) gives a Lagrangian writing of the energy balance in the local form,

d s dt ∫ Ω t E s + 1 2 E s f dΩ t + d f dt ∫ Ω t E f + 1 2 E sf dΩ t = ∫ Ω t d s E dt + E∇ • s + ∇ • 1 ρ f n E f + 1 2 E sf w dΩ t , = ∫ Ω 0 s dE dt + ∇ s • 1 ρ f n E f + 1 2 E s f M dΩ 0 s , ( .
dE dt = S . . E + b f • M m f -∇ s • s p f M m f + ∇ s • s F -1 s • Σ d f T • F s M m f -∇ s • ∇ s • γ f M m f - γ f s ∇ s • s M m f -∇ s • 1 ρ f n E f + 1 2 E sf M + q .
( . )

The second principle of thermodynamics e second principle of thermodynamics entails a fundamental inequality related to entropy evolution thus establishing a distinction between reversible and irreversible processes. For a porous medium saturated with a non-uniform uid this inequality can be wri en as,

d s dt ∫ Ω t ρ s (1 -n)s s dΩ t + d f dt ∫ Ω t ρ f ns f dΩ t ≥ - ∫ ∂Ω t q • n T dS t , ( . )
where s s has been introduced as the speci c entropy of the solid matrix.

Employing similar steps as used for the rst principle, the pull-back into the solid reference con guration of the entropy inequality reads in its local form,

dS dt + ∇ s • s f M + q T ≥ 0, ( . )
where S = s ρ s (1 -n)s s + ρ f ns f denotes the Lagrangian overall entropy per unit volume of the porous medium.

The Clausius-Duhem inequality & dissipation

Introducing the overall Helmholtz free energy of the porous medium, Ψ, though the Legendre transformation ) one can transform the rate of overall entropy within Eq.( .) to obtain,

Ψ = E -T S, ( .
dE dt -S dT dt - dΨ dt + T ∇ s • s f M + q T ≥ 0. ( . )
Using the expression for rate of overall internal energy from Eq.( . ) in Eq.( .) once can obtain the Clausius-Duhem inequality,

S . . E + b f • M m f -∇ s • s p f M m f + ∇ s • s F -1 s • Σ d f T • F s M m f -∇ s • ∇ s • γ f M m f - γ f s ∇ s • s M m f -∇ s • 1 ρ f n E f + 1 2 E sf M + q + T ∇ s • s f M + q T -S dT dt - dΨ dt ≥ 0.
( . )

Following a classical approach owing to the additive nature of energy and entropy, the Lagrangian densities of free energy, Ψ s , and entropy, S s , of the porous solid are deduced by subtracting respective bulk contributions of the uid from those of the overall porous medium,

Ψ s = Ψ -m f ψ f , S s = S -m f s f . ( . )
As mentioned at the beginning of this chapter, the above terminology is due to the work of B ( ), according to which the porous solid is understood as a 'we ed' porous skeleton, with a thin layer of immobile uid a ached to the pore walls and thus accounting for the associated solid-uid interfaces.

is allows one to model Ψ s as a state function in the context of thermoporoelasticity, where dissipation is only a ributed to ow of mobile uid and thermal e ects.

Introducing Eq.( . ) into Eq.( . ), while employing the balance of uid mass, Eq.( . ) and constitutive prescription of the non-uniform uid, Eq.( . ) the Clausius-Duhem inequality can be rephrased identifying an the overall dissipation of the porous medium Φ as,

Φ = - 1 S r ∇ s P + ∇ s ∇ s • γ f ϕS r + 1 S r 1 2n E s f + γ f ϕ • ∇ s s s - b f ϕS r M ρ w + S . . E + P + 1 2n E s f + γ f ϕ • ∇ s s s dϕ dt + ϕ S r 1 2n E sf + γ f ϕ • ∇ s s s dS r dt - γ f ϕS r • d∇ s (ϕS r ) dt -S s dT dt - dΨ s dt - q T ∇ s T ≥ 0. ( . )
As a general practice in poromechanics (C , ) the individual contributions by the solid skeleton (Φ s ), non-uniform uid ow (Φ f ) and thermal e ects (Φ T ), to the overall dissipation are identi ed and are independently assumed non-negative in order to satisfy the above inequality so that,

Φ s = S . . E + P dϕ dt + 1 S r 1 2n E s f + γ f ϕ • ∇ s s s d(ϕS r ) dt ( . ) - γ f ϕS r • d∇ s (ϕS r ) dt -S s dT dt - dΨ s dt ≥ 0, Φ f = - 1 S r ∇ s P + ∇ s ∇ s • γ f ϕS r + 1 S r 1 2n E sf + γ f ϕ • ∇ s s s - b f ϕS r M ρ w ≥ 0, ( .
)

Φ T = - q T ∇ s T ≥ 0. ( . )
Eq.( .) in addition to the postulate of local state leads to the free energy density of the solid being a rate independent function of the state variables: solid strain, Lagrangian porosity, saturation degree, spatial gradient of Lagrangian water content and absolute temperature, so that we obtain,

Ψ s Ψ s (E, ϕ, S r , ∇ s (ϕS r ),T ).
( . )

e dissipation due to uid ow is encompassed within Eq.( .). e essence of this restriction is that the generalized force, coe cient of M/ρ w , causing the uid ow must be related to the uid velocity vector in such a way as to satisfy the dissipation inequality. Eq.( . ) as usual states that heat ows from regions of higher temperature to those of lower temperature along the direction of negative spatial gradient of T .

. . Poroelastic constitutive relations

In what follows we assume isothermal conditions. Consequently the thermal dissipation is vanishing and the dependency of free energy density of the solid in Eq.( . ) reduces to,

Ψ s Ψ s (E, ϕ, S r , ∇ s (ϕS r )). ( . )

Solid state equations

Poroelasticity amounts to vanishing solid dissipation and absence of any frozen contributions to the free energy. Within this framework, using Eq.( . ) in Eq.( .) gives the following restriction on Ψ s ,

Φ s = S - ∂Ψ s ∂E dE dt + P + 1 2n E s f + γ f ϕ • ∇ s s s - ∂Ψ s ∂ϕ dϕ dt + ϕ S r 1 2n E s f + γ f ϕ • ∇ s s s - ∂Ψ s ∂S r dS r dt - γ f ϕS r + ∂Ψ s ∂∇ s (ϕS r ) • d∇ s (ϕS r ) dt ≥ 0, ( . )
wherein variations of the state variables can occur independent of each other. Similar to classical approach, the following set of dual relations can be obtained as state equations of the porous solid,

∂Ψ s ∂E = S, ( .
)

∂Ψ s ∂ϕ = P + 1 2n E s f + γ f ϕ • ∇ s s s , ( . ) ∂Ψ s ∂S r = ϕ S r 1 2n E s f + γ f ϕ • ∇ s s s , ( .
)

∂Ψ s ∂∇ s (ϕS r ) = - γ f ϕS r . ( . )
Eq.( .) is classical in poromechanics giving rise to a constitutive relation for total stress tensor. Eqs.( . , . ) are reminiscent of the standard relations in unsaturated poromechanics (see Sec.( . . ) of C ( )) with additional contributions due to the hyper stress of the uid, which itself is restricted through Eq.( . ).

In fact, in the absence of gradients of uid content, Eq.( .) reduces, allowing a direct comparison with the classical relation involving macroscopic capillary pressure, p c , (Eq.( .)) as follows,

∂Ψ s ∂S r ≈ ϕ S r 1 2n E s f = -ϕp c = ϕ ∂U ∂S r . ( . )
Assuming such comparison to be valid, Eqs.( . , . ) can be rephrased introducing p c ,

∂Ψ s ∂ϕ = P + -S r p c + γ f ϕ • ∇ s s s , ( .
)

∂Ψ s ∂S r = -ϕp c + γ f ϕ • ∇ s s s . ( . )

Generalized Darcy's law

Coming to the uid dissipation, Eq.( . ), a way to satisfy this restriction, which is a usual practice in poromechanics, is to assume that the dissipation is a quadratic function of uid velocity vector, V = M/ρ w , the coe cient of proportionality being the inverse of the second-order permeability tensor,

K. With such an assumption, M is constitutively restricted to

M = -ρ w K • 1 S r ∇ s P + ∇ s ∇ s • γ f ϕS r + 1 S r 1 2n E sf + γ f ϕ • ∇ s s s - b f ϕS r . ( . )
Introducing Eqs.( . , . , . & . ) into the above results in a generalization of the unsaturated Darcy's law,

M = ρ w V = -ρ w K • ∇ s ∂Ψ f ∂S r + ∂U ∂S r -∇ s • ∂Ψ N L ∂∇ s (ϕS r ) - b f ϕS r .
( . )

In the case of isotropic porous medium the permeability tensor reduces:

K = (κ/η w )K(S r )I , where K(S r
) is called a ux function and is assumed to play an equivalent role as the relative permeability of we ing uid k w (S w ) in classical partial saturation introduced in Sec.( . . ). is equivalency in part is justi ed due to the choice of the phase eld parameter S r as the saturation degree of the we ing uid. With the above setup, the uid dissipation in Eq.( . ) now reads,

Φ f = η w κ K(S r ) (M • M) ρ 2 w ≥ 0. ( . )
is results in a restriction on the function K(S r ) to be positive in order to ensure the positiveness of the dissipation associated to the uid ow. Moreover, it is worth to note at this point that there is no restriction due to thermodynamic principles on the convexity of the function K(S r ). e particular form employed for this function in the current work is motivated further in Chp.( ).

Towards a non-uniform pore fluid

Looking at Eq.( . ), in the absence of a bulk force acting on the uid, b f = 0, the generalized force driving the ow is the negative spatial gradient of a chemical potential identi ed as,

µ pf = ∂Ψ f ∂S r + ∂U ∂S r -∇ s • ∂Ψ N L ∂(∇ s (ϕS r )) , ( .
)
where a new sub-script 'p f ' is introduced that represents a supposed 'non- uniform pore uid' that is understood as a non-uniform uid, as introduced in Sec.( . . ), that is con ned within a porous network characterized by the capillary energy U (S r ). As usual in Cahn-Hilliard type modeling, this chemical potential, µ pf , can be derived from the variational derivative w.r.t S r of an overall free energy of the non-uniform pore uid when ϕ = const, i.e., when the porous skeleton is rigid. e density of such an overall free energy, Ψ pf , in the current case can clearly be identi ed by the overall free energy density of the porous medium, Ψ, assuming a rigid porous skeleton given by,

Ψ pf = ϕ Ψ f (S r ) + ϕ U (S r ) + Ψ N L (∇(ϕS r )). ( . )
It is worth noting that unlike Ψ f (S r ), the local part of this new energy density, (Ψ f (S r ) + U (S r )), has no more a symmetric double-well structure between S r = 0 and 1. Instead, the minimum associated to S r = 0 no more exists and only one global minimum remains at S r = 1. See Fig. ( . ). Moreover, depending on the relative intensities of bulk and capillary energy densities it is possible also that a minimum associated to lower S r exists, but is shi ed inwards of the original range of S r thus changing the corresponding local preferential states at equilibrium. However, these two minima would not be isopotential. In both these cases, when the chemical potential is non-monotonic, Maxwell construction (equal-area rule) can be understood as a search for the line bi-tangent to this new energy density at equilibrium conditions of the non-uniform pore uid. is is equivalent to construction of a linear potential due to external forces which once accounted for brings back the double-well structure with two isopotential minima. Now these two minima, both shi ed inwards of original range of S r , being isopotential allows for coexistence between the corresponding phases. e inward shi of the minima can be understood as a correction to account for we ing properties of the skeleton at equilibrium in the case of lower minimum and as a correction to account for trapped air in the case of higher minimum shi ed from S r = 1. Owing to this signi cance of physical interpretation, this local contribution in Eq.( .) is further referred to as the e ective energy density, Ψ e (S r ) = Ψ f (S r ) + U (S r ), of the non-uniform pore uid and its partial derivative w.r.t S r , the e ective chemical potential, see Fig. ( . ),

µ e (S r ) = ∂Ψ f ∂S r + ∂U ∂S r = µ f (S r ) -p c (S r ). ( . )
. .

Extended Richards' equation: Rigid porous skeleton

Now that a framework resulting in thermodynamic restrictions has been recalled in su cient detail, we move forward with simplifying assumptions.

Restricting ourselves to the scope of the current thesis, we assume a rigid porous skeleton (ϕ = ϕ 0 = const) that is saturated with the non-uniform pore uid. Such an assumption encompasses the hypothesis of small perturbations as already introduced in Sec.( . . ). Moreover we assume isotropic linear properties of the porous medium.

Assumption of a rigid porous skeleton allows one to neglect the power of external forces acting to deform the skeleton itself and consequently any power of internal forces acting to generate solid strains. So, any dissipative phenomena should be restricted to the uid ow and interfacial changes. It is to be noted that with the introduction of the non-uniform pore uid, the interfacial energy contributions, which were previously accounted for within the free energy of the porous solid, Eq.( . ), are now part of Ψ pf , given by Eq.( . ). is allows to describe the dissipation within the system solely by ow of this non-uniform pore uid through Eq.( . ) and the evolution of S r using the mass balance, Eq.( . ). Both of these read respectively as follows under the said assumptions,

M = ρ w V = -ρ w κ η w K(S r ) • ∇ ∂Ψ f ∂S r + ∂U ∂S r -∇ • ∂Ψ N L ∂∇(ϕS r ) - b f ϕS r , ( . ) ϕ ∂S r ∂t + ∇ • V = 0. ( . )
Similar to the developments in Sec.( . . ), an extension to Richard's equation, Eq.( . ), can be derived in the current phase eld context by substituting Eq.( . ) into Eq.( . ),

ϕ ∂S r ∂t -∇ • κ η w K(S r ) • ∇ ∂Ψ f ∂S r + ∂U ∂S r -∇ • ∂Ψ N L ∂∇(ϕS r ) - b f ϕS r = 0.
( . )

e appropriate boundary and initial conditions that supplement the above equation to form the BIVP are detailed further in Part II in speci c contexts of imbibition and drainage.

. e current chapter serves as a basis for further developments in this thesis. e developments of Sec.( . ) are utilized for the study mechanical instabilities through damage gradient modeling in Part III. e extended Richards' equation developed in the context of the phase eld approach to partial saturation in Sec. ( . ) In the previous chapter a phase eld approach has been introduced to treat partially saturation and an extended Richards' equation has been derived in the context of non-deformable porous skeleton. A one-dimensional analysis and resolution of solutions is presented in the current chapter for the extended Richards' equation. One-dimensional PDE simulations reveal that the model is able to describe both imbibition and drainage fronts when employing a non-convex ux function. Further these moving fronts are resolved using a traveling wave analysis thus se ing the stage for a linear stability analysis against imposed perturbations which is the intent of the next chapter.

. e current analysis starts from the extended Richards' equation, Eq.( . ), derived under the assumption of a rigid porous skeleton (ϕ = ϕ 0 = const). In soil hydrology and various other practical applications like CO 2 sequestration, the bulk force acting on the uid is due to gravity. If we consider gravitational acceleration to be acting in the positive x-direction, then b f = m f e x ≈ ρ w S r ϕ e x in accordance with Eq.( .).

e structure of Eq.( . ) closely resembles that of the classical Richards equation, Eq.( .).

e di erence lies in the constitutive prescription that lead to the particular expression of µ pf , Eq.( . ), compared to that of p c (S w ) in Eq. ( . ).

In what follows, a more general case is considered by introducing an additional linear pressure distribution which potentially allows to describe an initial uniform background mean ow. is state can then be perturbed by modifying the boundary conditions in order to investigate the evolution of either a drainage or an imbibition front. us an augmented chemical potential incorporating both the e ects of initial background mean ow and of gravity forces is wri en as, µ = µ pf + P, ( . )

where the pressure distribution P = -λx. In the case when only gravity forces are considered P is just proportional to λ = ρ w . Both these scenarios may coexist in which case gravity can act to either stabilize or destabilize the motion of an air-water front depending on their relative directions.

Introducing the above mentioned augmented chemical potential, µ, the governing equation for S r , Eq.( . ), can be re-wri en as follows,

ϕ ∂S r ∂t + κλ η w ∂K(S r ) ∂x - κ η w ∇ • K(S r ) ∇µ pf = 0, ( . )
with µ pf given by Eq.( .) and reported again below for the sake of clarity,

µ pf = ∂Ψ f ∂S r + ∂U ∂S r -∇ • ∂Ψ N L ∂(∇(ϕS r )) . ( . )
e quantity κλ/η w [LT -1 ] can be identi ed as the magnitude of a saturated mean velocity, V m , in the porous medium. In essence this quantity represents the intensity of advection of non-uniform pore uid, in other words that of the phase eld parameter S r . Appropriate boundary conditions will be introduced in Secs.( . . & . . ) to close the problem.

. & e e ective dimension of Eq.( .) is [T -1 ]. To render the equation dimensionless the following dimensionless variables (.) and corresponding characteristic scaling numbers (sub-scripted 'h') are introduced:

x = x x h , μpf = µ pf p h = µ pf R Cγ , t = t t h = tV m x h . ( . )
e time scale is chosen to make unity the coe cient in front of advection term, t h = x h /V m . And a particular grouping of parameters results in the dimensionless equation,

ϕ ∂S r ∂ t + ∂K(S r ) ∂x - δ R C a ∇ • K(S r ) ∇μ pf = 0. ( . )
e dimensionless e ective chemical potential of the non-uniform pore uid, μpf , along with the linear pressure contribution, P, can be wri en, employing the particular expressions introduced in Secs.( . . & . . ), as,

μ = μpf + P = 2S r 1 -3S r + 2S 2 r -D U S r -S r es r 1 -S r es r -1 m -1 1-m -D N L ∇ • ∇(ϕS r ) - C a δ R x .
( . )

e dimensionless numbers C a , δ R , D U and D N L are identi ed as follows:

C a = V m η w R 2 κCγ , δ R = R x h , D U = π 0 R Cγ , D N L = C k R Cγ x 2 h = π k R Cγ x h 2 . ( . )
Here, C a is the Capillary number signifying the competition between viscous forces and capillary forces at the air-water interface. δ R is the ratio of characteristic radius of porous channels to the reference length scale. D U is a measure of relative strength of the coe cients of capillary interfacial energy and the air-water interfacial energy. And D N L is the so-called Cahn number, signifying intensity of the gradient energy with respect to the di usive term. Above in Eq.( . ) it is shown that the expression of D N L can be recast such that the interface thickness, , is introduced into the equations and a free parameter, π k [ML -1 T -2 ], appears. is allows an input to the model from experimental measurements of observed macroscopic transition lengths. Further in this work, the length scale, x h , is chosen to be the characteristic physical length of the domain under consideration. ese choices are shown further to produce transition lengths of order in the numerical solutions when ow is driven by gravity. e full dimensionless form using Eq.( . ) in Eq.( .), with the above choices of dimensionless numbers is as follows,

ϕ ∂S r ∂t + ∂K(S r ) ∂x - δ R C a ∇ • K(S r ) ∇ 2S r 1 -3S r + 2S 2 r -D U S r -S r es r 1 -S r es r -1 m -1 1-m -D N L ∇∆(ϕS r ) = 0.
( . )

It is to be noted that (.) has been dropped in Eq.( .) and is done so further in this work. So, from this point all the variables are dimensionless unless either mentioned otherwise or referred to from earlier sections.

.

&

In the work of S T ( ), a fundamental study of the stability of uid-uid interface has been done in an analogous Hele-Shaw ow context. In order to perform this, an horizontal sharp interface has been assumed to separate two uids of di erent viscosities moving within a Hele-Shaw cell, under the forces of gravity and pressure gradient. Upon this horizontal interface, wave like disturbances of variable wavelengths have been assumed and their growth in time has been understood as the typical ngering instability. Furthermore, surface tension has been shown to introduce a lower bound for the range of wavelengths of disturbances for which the interface is unstable.

It our intention here, to characterize within the framework of the adopted model, propagation of an air-water interface both in the case of imbibition and in drainage. And in the Chp.( ) the stability/instability of these interfaces will be investigated. erefore, as a rst step transversely homogeneous solutions of Eq.( .) that represent such interfaces need to be built.

In this section we observe that solutions assuming transverse homogeneity in and z-directions but evolving longitudinally along x-coordinate, are composed of similarity solutions of the one-dimensional equation

ϕ ∂S r ∂t + ∂K(S r ) ∂x - δ R C a ∂ ∂x K(S r )µ e (S r ) ∂S r ∂x + δ R C a ϕD N L ∂ ∂x K(S r ) ∂ 3 S r ∂x 3 = 0.
( . )

In the above, the de nition of e ective chemical potential of non-uniform pore uid, µ e (S r ) (Eq.( . )), is invoked. (.) denotes a partial derivative w.r.t S r throughout this Part II.

At larger spatial scales compared to both the length scale of the airwater interface and that of di usion, the solutions of Eq.( . ) are regularized solutions of the corresponding scalar hyperbolic conservation law (L F , ) in the limit of vanishing di usion. Assuming such limit, Eq.( . ) simpli es as

ϕ ∂S r ∂t + ∂K(S r ) ∂x = 0. ( . )
Here, the function K(S r ) plays the role of an advective ux. A class of weak solutions of a Riemann problem governed by Eq.( .), given a piece-wise uniform initial condition with a jump between S -and S + at x = 0, that represent sharp displacement of one uid by another, are piece-wise uniform functions known as shocks,

S r (x, t) = S -if x < ct S + if x > ct, ( . )
that move with a characteristic speed, c. ese solutions are self-similar in nature with respect to the transformation ξ = xct and propagate along the x-coordinate while satisfying the Rankine-Hugoniot jump condition,

K(S + ) -K(S -) = cϕ(S + -S -), ( . )
that relates the speed of the shock wave to the uniform solution values on either side of the shock. A shock is considered classical or 'compressive' if the characteristics on either side of the shock impinge onto it. is condition is given by the celebrated Lax Entropy condition, which in the current case can be wri en as

K (S + ) ϕc K (S -). ( . )
Apart from shocks, the class of smooth monotone weak solutions of the Riemann problem governed by Eq.( . ) are rarefaction waves connecting S -to S + . ese solutions are expansive and self-similar with respect to the transformation ξ = x/t. ese are described by the following form:

S r (x, t) =          S - if ϕx < t K (S -) K -1 (ϕx/t) if t K (S -) < ϕx < t K (S + ) S + if ϕx > t K (S + ).
( . )

When the ux function, K(S r ), is purely convex and increasing, the only possible weak solutions are either classical shocks (if S -> S + ) or rarefaction waves (if S -< S + ) and analogously for purely concave ux function (L F , ). In the presence of di usion these sharp classical shock solutions satisfying Eq.( . ), tend to be smeared up to a nite distance. On the other hand, for ux functions that have in ection points the structure of solutions is much more rich in the presence of higher order di usion/dispersion. For instance smeared shock solutions that violate Eq.( . ) can exist, which are considered to be non-classical (B ., ; H S , ). e class of non-classical shocks for which the characteristics on either side pass through the shock are termed 'undercompressive' (D ., ) and they satisfy either of

K (S ± ) ϕc ; K (S ± ) ϕc. ( . )
And those for which the characteristics behind and in the front seem to expand the shock itself are termed 'expansion' shocks (E ., ) and they violate Eq.( . ) as

K (S -) ϕc K (S + ).
( . )

Presence of these solutions and their manifestation as traveling waves in the solution structure of imbibition and drainage problems governed by Eq.( .) are shown in the following section. Further in Sec.( . ), the traveling wave part of these solutions are resolved in a one-dimensional se ing.

As mentioned earlier, the nature of the ux function, K(S r ), has an important e ect on the structure of solutions that represent air-water displacements. In classical partial saturation, Sec.( . ), the role of modulating the individual uid velocities in space accordingly to their respective saturation degree is played by the relative permeability function. In the current study we intend to use the functional form of such a relative permeability of water for the ux function, owing to the phase eld parameter, S r , of the current model being representative of the saturation degree of water. Typically unsaturated relative permeability functions used for various soils and uid combinations are parametric models in order to t experimental data. Especially in soil hydrology when uid combination is water and wet air the well known van Genuchten model is widely use as introduced in Sec.( . . ).

e expression of relative permeability of we ing uid is recapitulated below with a change of notations for convenience,

K(S r ) = S r 1 -(1 -S r a ) b c , (

.

In the current section we characterize the spatio-temporal evolution of solutions both in imbibition and drainage scenarios. To do this we choose the primary unknowns as the Saturation degree, S r , and the regularized e ective chemical potential with known spatially linear pressure contribution, µ. en the coupled system of equations formed by Eq.( . ) and Eq.( .) is resolved. We acknowledge at this point that since our focus is to analyze the general structure and evolution of the solutions, we have adopted a simplest numerical discretization. One can de nitely extend this to more sophisticated techniques of the likes of adaptive re nement (M ., ; B

., ) and non-local operator methods (R ., ).

Table . : Material properties, model parameters used through Part II, unless mentioned otherwise.

property/ parameter

Sand

κ [m 2 ] . E- η w [Pa.s] . E- ϕ [-] . C [-] . γ [N.m -1 ] . R [m] . E- π 0 [Pa] . m [-] . S r es r [-] . C k [N] .
[m] .

. E- δ R . E- D U . D N L . E-
Spatial discretization is done employing a standard Galerkin formulation and linear Lagrange nite elements. Time discretization is done using the implicit Euler scheme of rst-order. e discrete solutions at n th time step, S n r and µ n , are obtained by searching in the Hilbert space of admissible functions given by the Cartesian product H s × H µ , with

H s := S n r ∈ H 1 (Ω) : S n r = Sn r on ∂Ω s H µ := µ n ∈ H 1 (Ω) : µ n = μn on ∂Ω µ ( . )
de ned over the discretized domain Ω, so that the ordered couple (S n r , µ n ) solve the non-linear coupled variational system,

∫ Ω q ϕ S n r -S n-1 r ∆t dΩ + ∫ Ω δ R C a K(S n r ) (∇ q • ∇µ n ) dΩ - ∫ Ω N µ q δ R C a K(S n r )∇µ n • n dΩ N µ = 0, ∫ Ω S µ n -µ e (S n r ) + C a δ R x dΩ - ∫ Ω ϕD N L ∇ S • ∇S n r dΩ + ∫ Ω N s S ϕD N L ∇S n r • n dΩ N s = 0.
( . )

Here q and S are test functions belonging to the Hilbert space (H s ) 0 × H µ 0 of functions which vanish on ∂Ω s and ∂Ω µ where the values of S r and µ are speci ed respectively. n is the outward unit normal vector to the boundary where it is referred to. Ω N µ is the part of the boundary where normal derivative of µ is speci ed, which translates to imposing at that part of the boundary an injection or extraction velocity of the uid with a natural form,

V f Ω N s = - δ R C a K(S r )∇µ. ( . )
Ω N s is part of the boundary where normal derivative of S r is speci ed. For one-dimensional simulations the computational domain is chosen along the positive x-direction such that x ∈ [0, 1]. e corresponding physical length, L = 100m, is chosen to be su ciently large such that the analysis is satisfactorily close to the limit of vanishing di usion mentioned earlier. e material properties of the porous medium and the parameters of the model chosen for the purpose of demonstration are listed in Table .( . ), which are in the range typical of silica sands saturated with air-water mixture. e corresponding dimensionless numbers de ned in Sec.( . ) are listed in Table .( . ).

Mesh convergence behavior has been tested, once chosen initial and boundary conditions corresponding to an imbibition and a drainage problem, which are elaborated in Secs.( . . & . . ), to validate the numerical solution.

In particular successively re ned discretization steps have been considered. e results of this analysis are reported in Appendix A .

. . Imbibition

Displacement of air by water representing imbibition can be understood as a solution which transitions from a higher degree of saturation to a lower .( . ), the appropriate initial condition for imbibition is set to S 0 r = 0.20 throughout the domain, which is close to the residual saturation, S r es r , and µ 0 = µ pf (S 0 r ) + P. For all t > 0, the normal derivative of S r is set to vanish at both the boundaries, x = 0 and x = 1. See schematic Fig. ( . ). At the le boundary an injection velocity, .). e di used interface in the solution connecting a higher value of S r to a lower value represents the in ltration of water into a uid-poor domain under the presence of gravity acting in the positive x-direction. e pro le within the domain, of the regularized e ective chemical potential, µ pf , follows that of the S r solution according to Eq.( . ), with a strong gradient e ect due to the interface. While all the saturation pro les are nonmonotonic in the vicinity of S 0 r , there exist overshoot and non-overshoot behaviors behind the invading front. e reasoning for presence or absence of such overshoots and the oscillatory behavior around S 0 r is reserved for analysis in the further sections.

V f (x =0) = - δ R C a K(S -)∇P = K(S -)e x , (
However, in all these solutions the transition from S -to S 0 r seems to translate self-similarly in space suggesting the presence of traveling wave type higher order approximations of shock solutions. is observation is justi ed by a transformation into a TW-coordinate, ξ = xct, where c is given by Eq.( . ) with S + = S 0 r , see Fig. ( . ). is results in the transition region of the solution to collapse into the vicinity of a single location, ξ = 0. e self-similar shock solution, Eq.( . ), of the hyperbolic equation, Eq.( . ), as well is plo ed in the same TW-coordinate in Fig. ( . ). In Sec.( . . ) these TW-solutions are resolved and are shown to be classical or compressive in the sense of shocks.

. . Drainage

Drainage of water by air is understood as the contrary of imbibition, which is a transition from lower saturation degree to higher, moving in the direction of the higher saturation. So the initial condition is chosen as S 0 r = 0.99 which is close to fully saturated condition, and µ 0 = µ pf (S 0 r ) + P. Owing to the phase eld parameter, S r , of the current model being representative of the saturation degree of water, an injection ux of the form Eq.( . ) at the boundary would not be appropriate to induce drainage. Instead the initial condition is perturbed at the le boundary, e boundary at x = 1 is drained with a Dirichlet boundary condition on µ, Eq.( . ), and the normal derivative of S r is set to vanish for all t > 0 similar to the imbibition case. ) and since S -< S + the rarefaction wave that is observed is expected. On the other hand, it can be observed that for the same transformation the part connecting S a to S 0 r does not collapse into a similarity type solution. In fact, the ux function is non-convex within the interval (S a , S 0 r ) = (0.847, 0.99) and so a rarefaction wave with S -= S a and S + = S 0 r is not an associated weak solution.

At the junction of the expanding part of the solution and the uniform state, S b , longitudinal oscillations are observed which seem to grow with time, see Fig. ( . ). It is explained in Sec.( . . ) with the help of dynamical systems analysis that these oscillations are triggered due to the nature of the equilibrium state associated to S b . e growth of such oscillations with time on the other hand, is justi ed in Chp.( ) by analyzing the linear stability of one-dimensional uniform saturation states against longitudinal perturbations.

e overall solution structure represents a drainage scenario in which ahead of the macroscopic air-water interface, the air phase can start invading the pores. And behind such interface there exists a uniform saturation state of the water phase corresponding to S b , that is not displaced. Now, for S d r = 0.20, we vary λ, see Fig. ( . ). For lower values of λ = 0.5ρ w , 1.5ρ w , and thus of C a , the solution structure is similar to that of λ = ρ w with two intermediate uniform states S a and S b appearing at later times. However, for relatively higher values of λ = 150ρ w , 1500ρ w a single uniform state, S b ≈ 0.797, 0.791 respectively, appears which connects to the right directly to S 0 r through a sharper transition representing the air-water front. is intermediate state, S b , connects to S d r to the le through an expanding part. e regularized e ective chemical potential, µ pf , follows a similar structure in accordance with Eq.( . ).

A transformation into the TW-coordinate ξ = xct, where c is given by Eq.( . ) with S -= S b and S + = S 0 r reveals the self-similar nature of the sharper transition, see .

In this section we intend to resolve the TW part of the solutions observed in Sec. ( . ). With this intent we make a transformation of the independent variables into a moving coordinate,

S r (x, t) = s(x -ct) = s(ξ ) ∀x ∈ ; ∀t > 0. ( . )
Where the TW-solution, s(ξ ), and speed, c, are a priori not known. e TW-coordinate, ξ , is assumed to have the same speed as that of a shock representative of the transition within the solution S r (x, t). And so the TW-solution is independent of time and the transition within it is always centered at ξ = 0. Introducing Eq.( . ) into Eq.( . ) yields an ordinary di erential equation (ODE) that s(ξ ) needs to satisfy:

-ϕc ds dξ + dK(s) dξ - δ R C a d dξ K(s)µ e (s) ds dξ + δ R C a ϕD N L d dξ K(s) d 3 s dξ 3 = 0, ( . )
e existence of such traveling wave type similarity solutions for Eq.( . ) can qualitatively explain part of the complete saturation solution, S r (x, t), ranging between the two uniform states (see gures in Sec.( . )). In other words the TW-solution, s(ξ ), of Eq.( . ) can represent a branch of the PDE solution S r (x, t), with the corresponding boundary conditions on an in nite domain:

s (ξ =+∞) = s + , ds dξ (ξ =+∞) = 0, s (ξ =-∞) = s -, ds dξ (ξ =-∞) = 0.
( . )

Integrating Eq.( . ) once w.r.t ξ , using the boundary conditions at ξ = +∞ and assuming higher derivatives of the solution vanish as ξ → +∞ yields in a canonical form:

d 3 s dξ 3 = C a δ R ϕD N L K(s) ϕc(s -s + ) -K(s) + K(s + ) + δ R C a K(s)µ e (s) ds dξ .
( . )

e Rankine-Hugoniot jump condition, Eq.( . ), for the speed, c, of the representative shock can be recovered by employing Eq.( . ) as ξ → -∞. Since we are looking for TW-solutions satisfying boundary conditions Eq.( . ), the states s -and s + are two known equilibria of Eq.( .). In order to analyze their properties and the possibility of existence of other such equilibria it is convenient to rewrite Eq.( . ) into a system of rst-order ODEs:

s ξ = , ξ = w, w ξ = C a (ϕc(s -s + ) -K(s) + K(s + )) δ R ϕD N L K(s) + µ e (s) ϕD N L .
( . )

In the above sub-script ξ is used to represent derivative of a variable w.r.t to ξ . In the sense of a dynamical system with ξ playing the role of a time-like independent variable, the TW-solutions can viewed as trajectories connecting associated equilibria, (s e , 0, 0), of such system. e local stability properties of such equilibria determine the nature of the connection between them, (B S ,

). e eigen values of a corresponding linearized system in the vicinity of an equilibrium reveals the behavior of a small perturbation to it. ese eigen values are the roots of the depressed cubic equation,

β 3 - µ e (s e ) ϕD N L β - C a (ϕc -K (s e )) δ R ϕD N L K(s e ) = 0. ( . )
A detailed deduction of the above linearization is done in Appendix A .

According to the classical Cardano-Tartaglia formula for roots of depressed cubic equations, for µ e (s e ) > 0, if the discriminant of Eq.( .) is negative, then all three eigen values corresponding to the equilibrium are real. is condition can be simpli ed to,

C a ϕD N L δ R (ϕc -K (s e )) K(s e ) < 4 27 (µ e (s e )) 3 1 2 . ( . )
Else, only one eigen value of the three is real and the other two are complex conjugates giving rise to oscillatory behavior in the vicinity of such equilibria. On the other hand, for a µ e (s e ) < 0, one eigen value of the three is real and other two are complex conjugates, irrespective of the sign of the discriminant. ese relations form the basis for understanding the behavior of TW-solutions in the vicinity of their respective equilibrium states, (s e , 0, 0).

. . Imbibition

As observed in Sec.( . . ), the imbibition solutions are self-similar in the whole of the domain such that the uniform saturation states observed in the PDE solutions correspond to the imposed boundary conditions, S -, on the le and initial conditions, S 0 r , on the right. Hence, saturation states at the boundaries of TW-solutions in Eq.( .) are chosen as s + = S 0 r and s -= S -. And then the third order ODE Eq.( .) is numerically resolved over a domain of nite physical length, L = 100m, that is su ciently large in comparison to the characteristic interface thickness, , such that we get a satisfactory approximation of an unbounded domain. Within this domain, the transition region from s -to s + would be centered at ξ = 0. We discretize this domain using a uniform nite di erence grid and use a second-order accurate central di erence scheme to approximate the derivatives in Eq.( .) and ( . ). en the solutions are obtained solving the non-linear problem starting from an initial guess that is a sharp transition at ξ = 0 between the intended saturation values on either side. Fig. ( . ) shows these solutions and their corresponding transitions as connections on the graph of ux function, K(S r ), for various values of s -and λ = ρ w . e slope of these connections is the speed, c, of the corresponding shock that is represented by the TW-solution. As can be inferred from the plot, the relation between characteristic speed on either side of the connection, K (s -), K (s + ), and the speed of the shock is such that these connections are representative of classical compressive shocks that satisfy the Lax entropy condition, Eq.( . ). Moreover, similar to observations made in Sec.( . . ), there exist non-monotonicities of overshoot type and oscillatory type. ese solutions are qualitatively comparable to the results of one-dimensional experiments done by D C ( ) for constant rate in ltration of water into initially dry sand. In those experiments, while for the highest imposed ux the water saturation pro le was monotonic, for intermediate uxes non-monotonicities were observed, which grew stronger as the ux is lowered within the intermediate range. It can observed in Fig. ( . ) the same trend as in the experiments, of the non-monotonicities with respect to varying ux, Eq.( . ).

In order to aid the analysis of the behavior in the vicinity of equilibria, (s -, 0, 0) and (s + , 0, 0), these TW-solutions understood as trajectories running between the equilibria are plo ed, see Fig. ( . ), in the phase space governed by the third order ODE system Eq.( .). For all s -the solution close to s + is oscillatory. However, behind the front, for certain values of s -= 0.85, 0.80, there is a non-oscillatory overshoot and for s -= 0.75, 0.70 the overshoot is oscillatory as it leaves the equilibrium at s -. For s -= 0.96 there is neither an overshoot nor oscillations behind the front.

e oscillatory behaviors at both the equilibria can be explained analyzing the corresponding linearized eigen values. For all le equilibrium states, (s -, 0, 0), such that µ e (s -) < 0 results in oscillations, when leaving this state as seen in the case of s -= 0.75, 0.70. e right equilibrium state, (s + = 0.20, 0, 0), is such that µ e (s + ) > 0. However, the condition on discriminant being negative, Eq.( . ), is not satis ed for the choice of C a and D N L resulting in oscillations when arriving s + = 0.20.

e overshoot behavior on the other hand can be explained by employing singular perturbation techniques (W , ). To do this we introduce a small nonzero parameter, ϵ ≈ O( /x h ), and expand the solution using a perturbation series, s(ξ ) = s 0 (ξ ) + ϵs 1 (ξ ) + ϵs 2 (ξ ) . . . , which separates Eq.( . ) into a cascade of problems each governing the solution at a particular order. As ϵ → 0 the leading order solution, s 0 (ξ ), is governed by

-ϕc ds 0 dξ + dK(s 0 ) dξ - δ R C a d dξ K(s 0 )µ e (s 0 ) ds 0 dξ = 0. ( . )
Integrating once w.r.t ξ and employing boundary conditions, Eq.( . ), at ξ = +∞ gives,

-ϕc(s 0 -s + ) + K(s 0 ) -K(s + ) - δ R C a K(s 0 )µ e (s 0 ) ds 0 dξ = 0. ( . )
is results in a smooth implicit solution valid outside the transition zone and so-called the 'outer' solution:

ξ (s 0 ) -ξ m = ∫ ξ ξ m dΞ = ∫ s 0 s m δ R C a K(S)µ e (S) -ϕc(S -s + ) + K(S) -K(s + ) dS, ( . )
where ξ m is a reference value, chosen to be 0.0, at which an arbitrary s m ∈ (s + , s -), chosen to be 0.5, occurs. is implicit solution is plo ed in Fig. ( . ) for values s -= 0.96, 0.85 and s + = 0.20. It is clear that in both cases the solution is non-monotonic in ξ for a range of s 0 within (s + , s -) and thus non-physical. Whereas, the parts of the solution for s 0 within (s -, 1.0) and (S r es r , s + ) increase monotonically with s 0 . Now, in order to resolve the nonphysicality a weak solution can be constructed by traversing the non-physical region with a discontinuity. However, this jump in the saturation degree needs to comply with a continuity in the pressure whose gradient should exist. is condition can be met in the current model thanks to the doublewell structure of the energy density which as explained earlier gives rise to isopotential saturation states. For the choice of parameters in the current study these isopotential states can be obtained as (S c 1 , S c 2 ) ≈ (0.91, 0.23) using Maxwell construction. Now, exploiting the translational invariance of the implicit solution, Eq.( . ), a weak implicit solution of the following form can be built using (S c 1 , S c 2 ):

ξ (s 0 ) -ξ m = ξ (s 0 ) -ξ (S c 1 ), ∀ ξ < ξ m , ξ (s 0 ) -ξ (S c 2 ), ∀ ξ > ξ m .
( . )

It is to be noted that the choice of boundary conditions on saturation degree, (s -, s + ), is independent of the isopotential saturation states determined (S c 1 , S c 2 ). us a combination of solutions is possible and accordingly the appropriate part of the outer solution, Eq.( . ), needs to be chosen to further translate and construct the weak outer solution, Eq.( . ). For instance, if s -> S c

1 and s + < S c 2 the part of outer solution, ξ (s 0 ), that is to be translated in the construction of ξ (s 0 ) is the non-monotonic part within (s + , s -). See for example the case of s -= 0.96 in Fig. ( . (a)). On the other hand if s -< S c 1 and s + < S c 2 , the monotonically increasing outer solution, ξ (s 0 ), within (s -, 1.0) is shi ed to the right giving ξ (s 0 ) -ξ (S c 1 ), and the non-monotonic part within (s + , s -) is shi ed to the right giving ξ (s 0 ) -ξ (S c

2 ). See the case of s -= 0.85 in Fig. ( . (b)). e jump part of the weak implicit solutions thus constructed traversing the non-physical zone can be seen as do ed vertical lines in Fig. ( . ). In this manner the presence of overshoots in the PDE solutions in Sec.( . . ) and in the TW-solutions can be rationalized comparing the le end boundary condition of saturation with the isopotential saturation state closest to it, S c 1 . Similarly two other combinations of solutions can be envisaged when s + > S c 2 , which involve undershoot when connecting to equilibrium state at s + .

On the other hand, the structure of the 'inner' transition region can be resolved by making a transformation into a stretched variable, ξ = (ξξ m )/ϵ, thus restricting the problem into the transition region. Under such transformation and further expanding the resulting solution, s, using a perturbation series w.r.t ϵ allows one to identify the leading order 'inner' problem,

-µ e (s 0 )

ds 0 d ξ + ϕ d 3 s0 d ξ 3 = 0, ( . )
for s0 . Integrating Eq.( . ) w.r.t ξ using the boundary conditions that allow appropriate matching of the inner and outer solutions on one side, i.e., s0 = S c 2 and derivatives of s0 vanish as ξ → ∞, gives the implicit inner solution,

ξ (s 0 ) -ξm = ∫ ξ ξm d Ξ = ∫ s0 sm 2 ϕ ∫ s0 S c 2 µ e ( S) -µ e (S c 2 ) d S -1 2 d S, ( . )
where ξm = 0 is a reference value where an arbitrary sm = 0.5 is chosen to occur. is inner solution approaches smoothly on the other side, ξ → -∞, to s0 = S c 1 , see Fig. ( . ). us the leading order inner solution is a ected neither by the form of ux function, K(s), nor by the boundary conditions chosen, Eq.( . ). Once the implicit solutions Eq.( . ) and Eq.( . ) are determined, a solution valid over the full domain can be constructed using asymptotic matching techniques with appropriate assumption on the region of overlap of those solutions. One of such possible assumptions involves a uniform matching,

s(ξ ) = s 0 (ξ ) + s0 (ξ ) -S c 1 , ∀ξ < ξ m , s 0 (ξ ) + s0 (ξ ) -S c 2 , ∀ξ > ξ m .
( . )

Figure . : Implicit inner solution, Eq.( . ), scaled back from the stretched variable space, ϵ ξ (s 0 ), shown as a solid line in the space (ξ (s 0 ), s 0 ). is solution varies between S c 1 ≈ 0.91 and S c 2 ≈ 0.23. e weak implicit outer solution, Eq.( . ), is shown for both the boundary conditions (s -, s + ) = (0.96, 0.20) (dashed) and (s -, s + ) = (0.85, 0.20) (dot-dashed). ese involve a jump between (0, S c 1 ) and (0, S c 2 ) which are shown as crosses. Solutions are shown in a restricted range of ξ for clarity. In the Eq.( . ) above, s 0 (ξ ) is the inverse mapping of ξ (s 0 ) in Eq.( .), with an abuse of notation replacing ξ with ξ noting that both take values from the real number line. s0 (ξ ) is the inverse mapping of the implicit inner solution that is scaled back from the stretched variable space, ϵ ξ (s 0 ). e accuracy of such solutions depends on the assumptions involved in the asymptotic matching procedure.

While, the weak implicit outer solution, Eq.( . ), indicates whether the overshoot and undershoot behaviors are to be expected or not in the numerical solutions, it is prudent to note that the integrand in Eq.( . ) diverges to -∞ while approaching both s -and s + . And this is irrespective of the relation between s -and s + , i.e., either s -> s + or s -< s + . us integrating backwards from s m ∈ (s + , s -) to min(s + , s -) always implies ξ (s 0 ) → ∞ and integrating forwards from s m ∈ (s + , s -) to max(s + , s -) always implies ξ (s 0 ) → -∞. So, this implicit outer solution, is only representative of imbibition solutions and as such drainage solutions are not possible in the context of the leading order outer problem, Eq.( . ). As we will see further, it is the second gradient term of the energy that is absent in the outer problem which is responsible for drainage solutions when K(S r ) is non-convex. show respectively the e ect of di erent le boundary conditions, s -, with a given initial condition, s + = 0.20 and of di erent initial conditions, s + , with a given le boundary condition, s -= 0.85. In the la er, for larger values of s + , the oscillatory non-monotonicity ahead of the interface are stronger and more spread-out along ξ , whereas the overshoot behind the interface tends to reduce slightly. Increasing λ tends to result in a reduced spread along ξ of the non-monotonicities within both behind and ahead of the interface. Also, the strength of these non-monotonicities tends to increase with λ.

. . Drainage

Contrary to the imbibition solutions in Sec.( . ), the drainage solutions are not self-similar throughout the domain. Taking cue from these observations, in what follows we consider one of the equilibrium states in the PDE solutions as a known. And then we use an iterative approach to nd the unknown equilibrium state, employing a shooting method to integrate forward or backward in ξ , the dynamical system Eq.( . ), until the target known equilibrium state is reached.

To do this in the case of lower values of λ we take as input from the D PDE results the higher equilibrium state, (S a , 0, 0), as the target and apply a small perturbation of the order E-to an initial guess of lower equilibrium state, (S 0 b , 0, 0), in the direction of the eigen vector corresponding to an unstable eigen value at (S 0 b , 0, 0). en we proceed to integrate forward Eq.( . ) using the MATLAB solver-ODE s (S R , ) which is based on a modi ed Rosenbrock method and features adaptive step size. In order to provide stopping criteria for the integration and direction for the iterations we observe two qualitatively di erent types of trajectories starting from an (S 0 b , 0, 0) in the vicinity of lower equilibrium state obtained in the D PDE results, (S b , 0, 0). See the trajectory reaches the target S a with a positive slope and then diverges towards full saturation. For higher values of S 0 b on the other hand, it reaches a local maximum missing the target and then moves to a local minimum before diverging to full saturation. e required le equilibrium state would form a trajectory that reaches the target with a zero slope and curvature, up to a numerical error. Assuming such trajectory exists for an S 0 b within the above two cases, we input them as starting points of iterations in a

Regula-Falsi method (G ,

). e stopping criteria for the iterations is chosen as a tolerance to be satis ed by the di erence between the target S a and the saturation value at which a zero slope achieved. is tolerance is set to E-in this work. For the case of higher λ, where only one intermediate uniform saturation state, S b , was observed in the D PDE results, we follow a similar approach and stopping criteria for the iterations, but with an a priori known target (S 0 r , 0, 0) = (0.99, 0, 0). Again we observe two types of trajectories starting in the vicinity of equilibrium state (S b , 0, 0). In both cases we observe an oscillatory behavior while approaching the target. See .). e connections for lower values of λ are representative of non-classical 'expansion' shocks that violate Lax entropy condition as Eq.( . ). Since these solutions are connections between two equilibria, corresponding local stability properties can be analyzed in accordance with Sec.( . ). In the current phase eld model these non-classical 'expansion' connections occur between a le equilibrium state, (s -, 0, 0), which has locally one unstable real eigen value, two stable complex conjugate eigen values and a right equilibrium state, (s + , 0, 0), which has locally two unstable real eigen values, one stable real eigen value. And hence there are no oscillations observed neither when leaving the le equilibrium state along the unstable eigen vector nor when approaching the right equilibrium state along the stable eigen vector. In Sec.( . . ) it is observed that longitudinal oscillations appear at the junction of a rarefaction wave and the uniform intermediate saturation state associated to S b . It can be inferred that these oscillations are due to the complex conjugate nature of the stable eigen pair in the vicinity of the le equilibrium state.

On the other hand, the connections for higher values of λ violate Lax entropy condition, Eq.( . ), while satisfying the rst relation in Eq.( . ) and so are representative of non-classical 'under-compressive' shocks which travel faster than the characteristic speed on either side of the shock. For such connections the system governed by Eq.( . ) has a third equilibrium, (s e3 , 0, 0), such that, ϕc(s e3 -s + ) -K(s e3 ) + K(s + ) = 0.

( . )

However the solution only represents a connection directly from (s -, 0, 0) to (s + , 0, 0) that have the same local stability properties with one unstable real eigen value, two stable complex conjugate eigen values. And the later are responsible for the small oscillatory non-monotonicity when the solution approaches (s + , 0, 0). See inset in Fig. ( . (d)). Similar 'under-compressive' shocks have also been observed as solutions of thin lm ow equations (B ., ; B S , ). In fact, the dynamical system formed by such equations has a similar structure to Eq.( . ). However, a signi cant di erence is that the counterpart of the coe cient of in Eq.( . ) 3 , µ e (s), in the case of thin lm ow equations is a constant instead of a function of the unknown eld that can change its sign as in the current context.

. Starting from the extended Richards' equation, we have presented an appropriate dimensional grouping and a one-dimensional analysis exploiting the hyperbolic structure of the governing conservation law. One-dimensional PDE simulations of the same have informed on the presence of similarity solutions for appropriate initial and boundary conditions. In particular the part of the solution which is a similarity solution of the traveling wave type was understood as the air-water displacement front. It has been demonstrated that both imbibition and drainage solutions are possible due to the non-convexity of the proposed ux function. Further with an intent to analyze their linear stability, the part of the solution that is a traveling wave type has been resolved by making a transformation into an appropriate moving coordinate system. ese solutions were analyzed by exploiting the underlying dynamical system. e solutions corresponding to imbibition driven by gravity have been found to be non-monotonic under certain boundary conditions and the reason for such non-monotonicities has been analyzed. In the current model, the double-well structure of the non-uniform pore uid energy has been found to cause overshoot type non-monotonicities depending on the boundary conditions chosen, whereas the oscillatory nonmonotonicities are due to the dynamical nature of the equilibria associated to the boundary conditions. e solutions of imbibition under the in uence of gravity have shown good qualitative agreement against one-dimensional inltration experiments, with the trend of the relationship between in ltration ux and saturation degree overshoot being preserved. Drainage solutions connecting unknown boundary conditions as well were resolved using an iterative shooting method starting from le boundary condition taken from the PDE solutions. In Chp.( ) we have done a one-dimensional analysis of the extended Richards' equation along with the resolution of traveling wave solutions. In the current study we present a linear stability analysis of its solutions which describe both imbibition and drainage. e essence of linear stability analysis performed is to understand if a small enough initial perturbation to the base solution grows or decays asymptotically in time by employing the principles of Lyapunov stability analysis. e results of stability against transverse perturbations sheds light on the sensitivity of the ow stability on injection ux, imposed pressure gradient and initial saturation degree. Further, a stability analysis against longitudinal perturbations reveals the instability of certain solutions that is associated to the non-monotonic structure of the uid bulk energy.

. S T ( ), have studied the stability of an horizontal sharp interface separating two viscous uids within a Hele-Shaw cell, driven by an imposed pressure gradient and under the in uence of gravity. It has been shown that wave like disturbances of arbitrary wavelengths can grow in time and eventually destabilize the interface if the motion of the interface is towards the more viscous uid and if the velocity is su ciently large.

Further the e ect of surface tension has been suggested, referring to the work of C . ( ), to be acting to limit the range of wavelengths of the disturbances to which the interface is unstable. A study of stability of the di used interface formed within the non-uniform pore uid, that is modeled using the Eq.( . ), is warranted, given the analogous nature of porous media ow with respect to Hele-Shaw ow. With this intention the base solutions, resolved as traveling waves are understood as plane wave solutions that are independent of the direction orthogonal to their propagation. Arbitrary wave like perturbations are then superposed to these base solutions and the growth in time of the former is understood as early time ngering phenomenon.

. .

Perturbed problem e base solutions in the physical domain are represented as S 0 (x, t) and they propagate along the x-direction. e physical solution itself is assumed to be a perturbed one, S(x, , t), that is composed of the base solution S 0 (x, t) at the leading order and superposed disturbances, resulting in a regular series expansion, S(x, , t) = S 0 (x, t) + ϵS 1 (x, , t) + ϵ 2 S 2 (x, , t) . . . , ( . ) where ϵ represents the magnitude of the disturbance. is form of the perturbed solution is inserted back into Eq.( . ) which is further linearized by ignoring quadratic and higher-order terms in ϵ. is leads to two problems, one of the leading order, O(1), and another perturbed problem of order O(ϵ).

is perturbed problem,

ϕ ∂S 1 ∂t + ∂(K 0 S 1 ) ∂x - δ R C a ∆ K 0 µ e0 S 1 + δ R C a ϕD N L ∇ • K 0 ∇∆S 1 + K 0 ∇∆S 0 S 1 = 0, ( . )
governs the evolution of a generic disturbance up to order O(ϵ). In the above and what follows, functions sub-scripted (.) 0 are to be understood as evaluated with the base solution. Noting that the base solution S 0 (x, t) is uniform along -direction and subsequently employing ξ = xct yields the O(ϵ) perturbed problem governing s 1 (ξ , , t) in the TW-coordinate,

ϕ ∂s 1 ∂t -ϕc ∂s 1 ∂ξ + ∂(K 0 s 1 ) ∂ξ - δ R C a s 1 d 2 dξ 2 K 0 µ e0 + 2 ∂s 1 ∂ξ d dξ K 0 µ e0 + K 0 µ e0 ∂ 2 s 1 ∂ξ 2 + ∂ 2 s 1 ∂ 2 + δ R C a ϕD N L s 1 d dξ K 0 d 3 s 0 dξ 3 + ∂s 1 ∂ξ K 0 d 3 s 0 dξ 3 + dK 0 dξ ∂ 3 s 1 ∂ξ 3 + ∂ 3 s 1 ∂ξ ∂ 2 + K 0 ∂ 4 s 1 ∂ξ 4 + 2∂ 4 s 1 ∂ξ 2 ∂ 2 + ∂ 4 s 1 ∂ 4 = 0.
( . )

Now since any generic disturbance can be resolved into its Fourier modes, modal stability analysis amounts to knowing the growth in time of individual Fourier modes composing the disturbance. If any one of these modes is found to be growing in time, the ow as a whole is to be understood as unstable since any physical disturbance would contain all of these modes. Following this approach a spatially periodic form is assumed for s 1 :

s 1 (ξ , , t) = e ik +σ t s(ξ ), ( .
)
where k is characteristic wave number of the disturbance in -direction, σ is the exponential growth factor in time (with an abuse of notation used for the Cauchy stress tensor) and s(ξ ) is the amplitude of the wave-like disturbance allowing variation in direction of propagation, ξ . Introducing Eq.( . ) in Eq.( . ) yields:

A d 4 s dξ 4 + B d 3 s dξ 3 + C d 2 s dξ 2 + D d s dξ + E s -σ s = 0, ( . )
a linear homogeneous ODE with spatially varying coe cients identi ed as,

A = - δ R C a D N L K 0 , B = - δ R C a D N L dK 0 dξ , C = - 1 ϕ δ R C a -K 0 µ e0 -2ϕD N L k 2 K 0 , D = - 1 ϕ K 0 -ϕc -2 δ R C a d dξ K 0 µ e0 + δ R C a ϕD N L K 0 d 3 s 0 dξ 3 -k 2 dK 0 dξ , E = - 1 ϕ - δ R C a d dξ K 0 µ e0 -k 2 K 0 µ e0 + δ R C a ϕD N L d dξ K 0 d 3 s 0 dξ 3 + k 4 K 0 + dK 0 dξ . ( . )
ese coe cients are dependent on the base solution, s 0 , its derivatives w.r.t ξ and the wave number of the disturbance, k. Since the base solution satis es boundary conditions s -and s + while approaching uniformly the boundaries ξ → -∞ and ξ → +∞ respectively, s needs to vanish uniformly at those boundaries such that the O(ϵ) perturbed solution is admissible. is results in the following boundary conditions that any s governed by Eq.( . ) needs to satisfy,

s (ξ =+∞) = 0 , d s dξ (ξ =+∞) = 0, s (ξ =-∞) = 0 , d s dξ (ξ =-∞) = 0. ( . )
. .

Discretization and Asymptotic growth/decay

In the current section for each base solution, the pair s and σ are determined numerically varying the value of k. We approximate using a second-order accurate central di erence scheme the derivatives in Eq.( .) and ( . ). e physical domain is the same as that is used to resolve the base solutions, with L 0 = 100m and corresponding dimensionless domain ξ ∈ (-0.5, 0.5).

Subsequently the discretized problem for s and σ can be posed as an eigen value problem, (º -σ ) s = 0, ( . ) where º is a sparse matrix encompassing the coe cients at respective nite di erence (FD) nodes and is an identity matrix of the same size. Since for a non-trivial s, |º -σ | = 0, the solution set for σ is given by the spectrum of º, represented σ º . Following modal stability analysis, if sup{ (σ º )} > 0, then the corresponding perturbation grows exponentially in time according to Eq.( . ). Whereas if sup{ (σ º )} < 0, then all the perturbations, s, decay exponentially with time. In order to numerically approximate this critical eigen value, σ s = sup{ (σ º )}, for a given wave number, k, a Krylov-Schur algorithm available in MATLAB (S , ) has been employed. is relation between k and σ s (k) is plo ed as a dispersion curve for each base solution, s 0 . Fig. ( . ) shows the dispersion relation for base solutions which represent imbibition for various values of s -, s + = 0.20 and λ = ρ w . e wave number of the fastest growing perturbation, k f , is the one associated to the peak of the dispersion curve. For a large value of s -= 0.96, representing large injection velocity according to Eq.( . ), σ s (k) is found to be negative for all values k analyzed. is indicates an exponential decay of the arbitrary perturbations applied to those base solutions and thus a stable imbibition front which is the associated plane wave. As s -is reduced the base solutions start to involve non-monotonicities and the growth rate σ s (k) is found to be positive for a range of k. is indicates that perturbations within those ranges can potentially grow in time, triggering ngering type instabilities. As s -is reduced further, s -0.36, the dispersion curves again indicate exponential stability of the base solutions.

is is in qualitative agreement with the experimental observations by S F ( ); D C ( ) and the related hypothesis that overshoot pro les of saturation result in an unstable we ing front in gravity driven in ltration (N ., ; E G , ; E ., ; D C , ). Also, indicated in Fig. ( . ) is the path followed by the wave number of the fastest growing perturbation, k f , as s -is reduced. e interpretation of this path is made clear in Fig. ( . ), with the understanding that the wave length associated to the fastest growing perturbation, λ(k f ) = 2π /k f , is representative of the spacing between the early time ngers. And the intensity of the non-monotonicity is re ected by (s max -s -), where s max is the peak value of the saturation degree a ained by the imbibition front. For intermediate range of injection velocities the intensity of overshoot in the base solutions is larger and the corresponding nger spacing is smaller, whereas the growth rate is larger. As the injection velocity either increases moving towards fully saturated conditions or decreases lower than the intermediate range, the intensity of overshoot reduces, the corresponding nger spacing increases and the growth rate decreases indicating a transition to the stable regime. is strongly indicates a correlation between the injection ux and the characteristics of the associated ngering instabilities as observed in experiments (G ., b; Y H , ). A similar observation has been done in C F J ( a). Also, it is interesting to compare the structure of the perturbation with maximum growth rate for a given wave number. Fig. ( . ) shows this normalized eigen function, s, determined by solving Eq.( . ) with k = 140 and the associated base solutions, s 0 , superposed. For all values of s -the eigen function reaches a maximum within the di used interface between s -and s + of the base solution. However, there exists a relative di erence in the peak values and the spread along ξ of these peaks among these eigen functions. As s -reduces the peak value reduces and the spread along ξ increases. is indicates a shi towards weaker and more spread-out perturbations. It can be inferred from this observation that for base solutions that consist a stronger non-monotonicity (lower s -), a relatively weaker perturbation is su cient to render the ow unstable and vice versa. imbibition for various values of s + , s -= 0.85 and λ = ρ w and the path followed by the wave number of the fastest growing perturbation, k f , as s + is increased. For high values of s + the base solution seems to be stable. e intermediate range of s + shows unstable behavior and as s + is reduced further the tendency seems to be towards stable behavior. However, it is to be noted that the lowest value of s + in this work is restricted by the residual saturation degree S r es r = 0.1567. is behavior is also plo ed in Fig. ( . ) showing the relation between the qualitative indicator of nger spacing, λ k (k f ), and s + . A qualitatively similar observation was made in point-source in ltration experiments by B . ( ), where water was injected into initially moist sand packs, albeit within a lower range of initial water content. It was observed that the ensuing nger sizes decreased slightly with increasing initial water content starting from a completely dry condition and a er a certain value, the nger sizes increased exponentially.

For k = 140 the normalized eigen function, s, is more spread-out along ξ and has a lower maximum value of the peak when s + is increased, see Fig. ( . ). is indicates that for larger s + while a weaker and more di used perturbation can invoke instability in the ow, it grows relatively slower in accordance with Fig. ( . ). (c)) that represent imbibition for varying values of λ = 0.5ρ w , ρ w , 1.5ρ w and with xed boundary conditions s + = 0.20 and s -= 0.85. For the range of λ studied it is clear that base solutions with a relatively narrower spread along ξ (larger λ) when leaving the equilibrium state, (s -, 0, 0) have a wider range of k for which the perturbations are unstable and also larger values of k f and associated peaks, corresponding to a relatively faster growing perturbations. To summarize the case of imbibition, linear stability analysis of the in ltrating front provides results consistent with experimental evidence. Moreover, it tells us that for a given imposed pressure gradient, stronger the perturbation needed to destabilize the front, faster would be its growth and consequently faster would be the destabilized nger. And this condition is achieved when the non-monotonicity of the base solution is less pronounced.

For the drainage case, ose solutions representative of non-classical 'under-compressive' shocks, which were observed in the higher range of λ, are found to be stable towards transverse perturbations. Whereas, the solutions representative of non-classical 'expansion' shocks, at the lower range of λ, are found to be unstable within a particular range of wave numbers, k, in spite of those base solutions being monotonic. Among the la er, for the range of λ studied, as λ is increased the range of unstable wave numbers increased along with the wave number of the fastest growing perturbation, k f . From Fig.( . ) for k = 150, it can be observed that for higher values of λ the eigen function has a lower peak and a wider spread along ξ representing a weaker perturbation that renders the ow unstable. It can be noted that this tendency is contrary to what was observed in the imbibition case.

While the subsequent growth of an initial disturbance is governed by non-linearities inherent to the model, in the context of linearized analysis it can be expected that the wave number of the ensuing ngering instability at late times would be close to that of the fastest growing perturbation, k f . Also the upper bound of the range of unstable wave numbers serves as a critical wave number beyond which any perturbations decay. In this sense, linear stability analysis serves as tool to identify the minimum size of the transverse dimension of a porous medium domain in order to be able to invoke and observe ngering phenomenon. Taking cue from the above results, two-dimensional PDE simulations in Chp.( ) were initialized within appropriately sized domains and with an imposed transversal perturbation superposing a few random wave numbers in the neighborhood of k f in order to accelerate the onset of instabilities. As can be observed from the structure of the base solutions of the TWproblem, there exist regions of uniform saturation degree ahead and behind the front. e solution within these regions is subject to instability in the presence of longitudinal perturbations such as the oscillatory non-monotonicities observed in base solutions in the vicinity of the front, see Fig. ( . ). Such instability triggers the formation of banks of uid perpendicular to the direction of motion of the front, which further destabilize transversely to form bubbles behind the front. Also, when a front destabilizes due to transverse perturbations, as analyzed in Sec.( . ), their ensuing growth forms ngers, see Chp.( ). Along the length of those ngers, there is also a region behind their tips within which the saturation degree is uniform. e solution within this region can as well destabilize due to longitudinal perturbations resulting in breakup of those ngers into droplets. is is akin to the now well-studied Plateau-Rayleigh instability that explains jets of liquid destabilizing and forming coalesced droplets (E V ,

). In order to analyze the onset of such instabilities in the current context, we consider a longitudinally perturbed solution, S(x, t), of Eq.( . ), in the vicinity of a constant and uniform base solution, S 0 , and expand it in a regular series, . ) about the base solution, with ϵ now being the magnitude of longitudinal perturbation. Introducing such perturbed solution into Eq.( . ), and further assuming that the perturbation is small compared to the base solution results in an O(ϵ) perturbed equation,

S(x, t) = S 0 + ϵS 1 (x, t) + ϵ 2 S 2 (x, t) . . . , (
ϕ ∂S 1 ∂t + K 0 ∂S 1 ∂x - δ R C a K 0 µ e0 ∂ 2 S 1 ∂x + δ R C a ϕD N L K 0 ∂ 4 S 1 ∂x 4 = 0, ( . )
that governs S 1 (x, t). Here as well, functions sub-scripted (.) 0 are to be understood as evaluated with the base solution. In order to analyze if a longitudinally harmonic perturbation would grow/decay asymptotically in time, the following form is assumed for the O(ϵ) perturbation,

S 1 (x, t) = e ikx +σ t , ( .
)
where k is characteristic wave number of the disturbance in x-direction and σ is its exponential growth factor in time. Introducing Eq.( . ) into Eq.( .) results in an expression for the growth factor in terms of the base solution and the wave number,

σ (S 0 , k) = 1 ϕ - δ R C a K 0 µ e0 k 2 + ϕD N L k 4 -ikK 0 . ( . )
Similar to what was noted in the Sec.( . ) for transverse perturbations, any longitudinal perturbation resulting in (σ ) > 0 is to be understood as growing exponentially in time leading to instability about the uniform base solution. Fig. ( . ), shows the graph of the relation Eq.( . ) in the space of (S 0 , k, (σ )). For certain values of S 0 there exist a range of wave numbers k for which (σ ) > 0. is region of instability, shaded grey in Fig. ( . (b)), is given by the relation,

µ e0 < -ϕD N L k 2 , ( . )
between the slope of the e ective chemical potential evaluated at S 0 and the square of the wave number. While the respective base solutions for both these cases exhibit oscillations when leaving the uniform equilibrium state (s -, 0, 0), it will be shown further using two dimensional simulations in Chp.( ) that only the oscillations in the case s -= 0.70 are within the range of k that trigger the instabilities behind the front.

.

In the current chapter a systematic linear stability analysis is presented of the solutions of the extended Richards' equation developed in Sec.( .). e solutions analyzed were both of imbibition and drainage. Transverse stability analysis of the imbibition solutions has indicated that they are sensitive to imposed in ltration ux as observed in the experiments. On the other hand for the drainage front a larger range of imposed pressure gradient was studied and the drainage front seemed to stabilize at very high pressure gradients. An analytical result was also presented for the stability of one-dimensional uniform base solutions against longitudinal perturbations in the direction of propagation of the front.

is result indicated on the possibility of destabilization of the solutions within a particular range of uniform values of the base solutions which could result in the formation of trapped regions of dissimilar saturation in higher dimensions. In the following chapter two-dimensional simulations are done for that verify the results of the above stability analysis. In the current chapter two-dimensional numerical simulations of the extended Richards' equation are presented with an intention to verify the results of stability analysis done in Chp.( ). e simulation results reveal uid ngering phenomenon for the boundary and initial conditions which are found to be unstable in the stability analysis against transverse perturbations. Further due to the double-well nature of the uid energy a rich structure of uid invasion is revealed as predicted in the stability analysis against longitudinal perturbations.

In this chapter we use numerical simulations as a means to examine the non-linear growth of perturbations to a base solution in a two-dimensional spatial se ing. To do this the two-dimensional domain is initiated with a eld of S r (x, ) involving a smooth jump along the x-direction and extending homogeneously in the -direction, that represents a transition from one uid phase to another. To this eld, perturbations dependent on -direction are superposed, that are sinusoidal in nature involving a few randomized wave lengths. is perturbed initial condition is representative of physically realistic conditions where the air-water interface in a porous medium is not independent of -direction owing to intrinsic heterogeneities. However, it is to be noted that these imposed perturbations are only spatial in nature and does not involve perturbations in the value of S r itself.

e numerical method and discretization scheme are the same as that of the one-dimensional PDE simulations in Sec.( . ) where the weak formulations are wri en in a general manner. e material properties and dimensionless numbers are as well utilized from Tables.( . & . ) respectively.

e details of the discretization domain, boundary and initial conditions for each case are detailed further.

. e computational domain used for imbibition simulations is a rectangle, such that x ∈ [0, 1] and ∈ [0, 0.3], whose physical dimensions correspond to a length of L = 20m and a height of H = 6m. Displacement of air by water representing imbibition can be understood as a solution which transitions from a higher degree of saturation to a lower one and moves in the direction of the lower saturation. In order to realize this, the domain is initialized with a solution S r (t = 0) involving a smooth jump along the x-direction, between V f .n = 0 (x,y=0.3)

V f .n = 0 (x,y=0)
the saturation degree representing the in ltrating water, S -, and that of an initial saturation degree of an almost dry domain, S 0 r , such that S -> S 0 r . At t = 0, the region of the domain initialized with S -is restricted to a thin layer at the in ltrating surface, (x = 0, ) that is perturbed transversely in the -direction with a few randomized wave lengths within the vicinity of the fastest growing wavelength, λ k (k f ), given by the linear stability analysis with boundary conditions s -= S -, s + = S 0 r . e initial condition for µ is such that µ(t = 0) = µ pf (S r (t = 0)) + P.

For all t > 0, the normal derivative of S r is set to vanish at all the boundary surfaces. At the le boundary surface, (x = 0, ), an injection velocity, is is a consequence of higher saturation degree within this fringe region compared to that of the initially dry region, S 0 r , see the ngers increases and their growth rate decreases, e ectively suppressing the ngering phenomenon. is is veri ed in the two-dimensional problem, see Fig. ( . ), where as S -in Eq.( .) is varied within the range [0.70, 0.96], the same trend is observed. For the largest injection velocity, S -= 0.96, the initial perturbations are completely suppressed and the imbibition front in ltrates the domain homogeneously. An important correlation is observed between these solutions and the results of the linear stability analysis in Sec.( . ). Speci cally, the number of ngers that eventually appear a er the non-linear growth for each S -, see Fig. ( . ), is strikingly close to the predictions of linear stability analysis which, in principle, is applicable only in the early stages of growth of the perturbations. If this correlation is universally true for solutions of models that exhibit non-monotonic dispersion curves is a question that could potentially have important consequences.

V f (x =0, ) = - δ R C a K(S -)∇P, (
For S -= 0.70 the region behind the in ltrating front which is expected to have a uniform saturation degree of 0.70 is perturbed longitudinally due to the oscillations as seen in the one-dimensional solution, see Fig. ( . ). ese oscillations are longitudinal in nature, as elaborated in Sec.( . ), and grow in time while interacting non-linearly with the in ltrating front. e result of this kind of instability and the non-linear interactions is the formation of bubble like trapped regions of lower saturation, see last row of Fig. ( . ).

. e computational domain used for drainage simulations is a rectangle, such that x ∈ [0, 1] and ∈ [0, 0.15], whose physical dimensions correspond to a length of L = 20m and a height of = 3m. Drainage of water by air is understood as the contrary of imbibition, which is a transition from lower saturation degree to moving in the direction of the higher So the initial condition, S r (t = 0), is chosen similar to that of imbibition, with a jump along x-direction in the degree S to S 0 r such that S < S 0 r . e initially saturated region of the domain is assumed to be close to full saturation, 0 r = 0.99. e region that is relatively S r ) = S -is restricted to a thin layer at the le boundary, (x = 0, ). Similar to imbibition, this initial condition is perturbed with randomly chosen lengths within the vicinity of the fastest growing wave length indicated by the linear stability analysis. e initial condition for µ is such that µ(t = 0) = µ pf (S r (t = 0)) + P. = µ e (S -) + P(x = 0, ).

( . ) e boundary at the right, (x = 1, ), is drained with a Dirichlet boundary condition on µ, Eq.( .), and the normal derivative of S r is set to vanish, similar to the imbibition case. ese set of boundary conditions have an e ect of inducing extraction of water at the right boundary that is consistent with the imposed pressure gradient, while the normal derivative of µ pf vanishes. e top, (x, = 0.15), and the bo om, (x, = 0), boundaries are impermeable and the normal derivative of S r is as well set to vanish. Fig. ( . ) shows the evolution of S r (x, ) for t > 0 of a perturbed initial condition, for the case of S -= 0.20, under the above boundary conditions and for λ = 1.5ρ w . e solution along the length of the domain is consistent with one-dimensional PDE solutions in Chp.( ), with a self-similar traveling wave part connecting a lower saturation degree, S b , to a higher, S a , which represents a non-classical 'expansion' shock and this is connected on either side up to the boundary by expanding parts of the solution. However, the traveling wave part of the solution along the length of the nger has a di erent lower value compared to that of within the gap between the ngers, see In fact, a one-dimensional traveling wave connecting S b to S a that correspond to the solution along the length of the nger would have a higher speed, approximated by Eq.( .), compared to that of within the gap between the ngers, which explains the higher speed of the advancing tip of the nger compared to its base.

As the solution evolves, formation and detachment of droplet like structures occurs between the ngers. is is due to the growth of longitudinal perturbations to the part of the solution, S r (x, = 0.0775) = S b ≈ 0.45, between the di used interface and the expanding region on the le . As noted in Chp.( ), the origin of longitudinal oscillations in this region is related to the nature of the equilibrium state associated to S b of the dynamical system formed by the ODE ( . ). ese oscillations tend to grow because S b ≈ 0.45 is clearly within the unstable region of uniform base solutions dictated by Eq.( .), see Sec.( . ). is is not true for the solution along the length of the nger because the corresponding S b ≈ 0.35.

Similar to the imbibition case, the e ective chemical potential, µ e , is nonmonotonic and the non-local chemical potential, µ N L , follows the convexity of S r , taking positive values just behind the di used interface as can be seen in Fig. ( . ). In order to verify the prediction of the linear stability analysis done in Sec.( . ), we present here the evolution of the solution S r for di erent values of the gradient of imposed pressure distribution, λ. Inline with the prediction, as λ is decreased from 1.5ρ w , the nger spacing is reduced resulting in thicker ngers that grow slower, see Fig. expanding rarefaction part of the solution behind the di used front does not seem to follow the same rate of decay resulting in persistence of the initial disturbance at the le end of this part of the solution.

. In this chapter two-dimensional numerical results were presented which were based on a nite element implementation. ese results verify those of the linear stability analysis with respect to nger sizes observed and stabilization of the front when expected.

e rich structure of the solutions is a ributed to the various energy contributions that encompass the non-uniform uid potential. While the results were inline with the existing experimental observations, the model predicts additional features that maybe observed during air-water displacements. Detailed experiments are warranted both in imbibition and drainage in order to investigate these predictions. Initiation and evolution of complex fracture networks is ubiquitous in geomaterials, especially during drainage. is chapter presents an application of gradient damage modeling in order to study initiation of fractures in the context of drainage induced partial saturation in porous media. e study starts with a modeling approach that envisages drainage induced mode-I fractures as caused by the build up of tensile stresses. Following this the evolution problem for damage is derived employing three principles underlying the variational approach to damage: irreversibility of damage, stability and energy balance. e governing equations of the problem are derived as the rst-order optimality conditions of the evolution problem. A numerical discretization scheme and algorithm are presented followed by a two-dimensional test case. e results portray the bifurcation and eventual localization of the damage into a series of well-spaced fractures originating at the drainage boundary and propagating into the drying porous skeleton driven by desaturation. e results also show qualitative similarity with observations in experimental desiccation tests. While the presented modeling approach remains a simplistic one, it is contended in Ch.( ) regarding its applicability to ne-grained soils which do not sustain signi cant tensile stresses when close to full saturation and an alternative underlying phenomenon is explored.

M E C HANICAL I N S TA BILI TIES

In this chapter the phase eld model of partial saturation is not adopted and the starting points of the current development are the classical poromechanical approach to partial saturation (see Sec.( . )), and the now well studied gradient damage modeling of bri le fracture (see M . (

) for an overview). In addition to this we work under an overarching assumption of a hydrophilic porous skeleton whose pore spaces are lled by two uids, air and water, between which air is considered passive and the hypothesis of small perturbations is assumed as detailed in Sec.( . . ).

.

Apart from the hydraulic problem related to the uid ow, the modeling developed in this chapter is done speci cally considering the porous skeleton. is is done adopting the hypothesis of energy separation as explained in Sec.( . . ). A free energy density of a partially saturated deformable porous skeleton is thus identi ed as a state function, G s (ε, π ). e corresponding state equations are the constitutive laws Eq.( . ) and Eq.( .), respectively for the dual variables total stress, σ (ε, π ), and Lagrangian porosity, ϕ(ε, π ), recapitulated below for easy reference.

σ -σ 0 = σ e -b(π -π 0 )I , with σ e (ε) = C . . (ε -ε 0 ), ( . ) ϕ -ϕ 0 = b(ϵ -ϵ 0 ) + 1 N (π -π 0 ), with π (p w , S w ) = p w S w -U (S w ). ( . )
Notwithstanding this modeling choice, it is to be noted that an equivalent formulation can be derived considering the porous solid (the so-called 'we ed' porous skeleton which comprises the skeleton and a thin layer of uid a ached to it) starting from Ψ s (ε, ϕ, S w ). is alternative approach is pursued in Ch.( ) to serve the investigation of a speci c phenomenon.

. . Gradient damage modeling setup

We consider an n d -dimensional domain, Ω ∈ n d , representing the initial con guration of a damaging isotropic partially saturated porous skeleton whose boundary is denoted by ∂Ω. Now, in line with standard gradient damage modeling (M ., ) and the theory of partially saturated poroelasticity, the following considerations are done:

. e scalar internal variable, α ∈ [0, 1], is assigned the role of describing the extent of damage. α = 0 denotes an intact healthy skeleton. Whereas, α = 1 denotes a fractured skeleton whose sti ness is degraded and can only sustain vanishing stresses.

. For a given unit volume of the porous skeleton in its reference con guration the bulk energy density, W s , is a state function characterized by: (ε, , ∇α), respectively the linearized strain tensor, the equivalent pore pressure, the damage variable and the gradient of damage.

. e skeleton is assumed to be dissipative in a non-local sense due to the dependency of the bulk energy density on the gradient of damage. e particular expression of the bulk energy density is assumed to be:

W s (ε, π , α, ∇α) = G s (ε, π , α) + w(α) + 1 2 w 1 2 ∇α • ∇α . ( .
e di erent terms of the above expression are explained below.

(a) With an abuse of notation G s (ε, π , α) now denotes the elastic energy density of the damaging porous skeleton which, as mentioned earlier, encompasses the contribution due to the deformation of the porous skeleton resulting in solid strains and changes in equivalent pore pressure accounting for the consequence of formation/annihilation of interfaces. A possible form this energy density that gives back the proper constitutive relations as its state equations reads as,

G s (ε, π , α) = 1 2 a(α)(ε -ε 0 ) . . C . . (ε -ε 0 ) -b(ϵ -ϵ 0 )∆π - 1 2N ∆π 2 -ϕ 0 π , ( . )
where the notations are borrowed from Sec.( . . ). ∆π represents the increment of equivalent pore pressure, (π -π 0 ). α 0 is a priori assumed to be a null eld associated to an initially undamaged skeleton. C is the sti ness tensor of the undamaged skeleton and its inverse, the compliance tensor, is denoted S. a(α) denotes a damage law that, when α grows from 0 to 1, degrades the elastic energy density responsible for solid strains to a vanishing value. While various choices of the functional form of a(α) are possible, we make a speci c following P . ( a),

a(α) = (1 -α) 2 + k , ( .
)
where k is a small positive constant used solely for numerical purposes to govern the fully damaged state. e consequence on the damage evolution due to the above choices is discussed further.

(b) e second term in Eq.( .) is the local part of the dissipated energy,

w(α) = w 1 α, ( . )
growing from 0 when α = 0, to a positive constant w 1 < +∞ when α = 1. us in accordance with the developments in P . ( a), since w (α) > 0 there exists an elastic phase preceding damage initiation and the niteness of w 1 ensures that the energy dissipated during a homogeneous evolution of α from 0 to 1 is as well nite.

e signi cance and possible motivation for the magnitude of this constant is discussed at a later moment.

(c) e last term in Eq.( .) is the non-local dissipation, which is assumed to be a quadratic function of the gradient of α, and is intended to regularize the local model allowing localizations of nite thickness. appears with the physical dimension of a length that in the context of gradient damage modeling is intended to control the localization thickness.

. e dual relations associated to the state variables in Eq.( .) are obtained as follows,

∂W s ∂ε = a(α)C . . (ε -ε 0 ) -b∆πI, ∂W s ∂π = -b(ϵ -ϵ 0 ) - 1 N ∆π -ϕ 0 , ∂W s ∂α = 1 2 a (α)(ε -ε 0 ) . . C . . (ε -ε 0 ) + w (α), ∂W s ∂(∇α) = w 1 2 ∇α . ( . )
According to the developments in P M ( ), the above choices of the functional forms of a(α) and w(α) result in a material behavior that is characterized by niteness of the energy dissipated during a full process of homogeneous damage and stress-so ening. is la er property of stress-so ening (resp. strain-hardening) requires that the elastic domain in the stress space (resp. strain space) is a decreasing (resp. increasing) function of α. ese spaces within the current context can be wri en employing the local part of the energy density, W s (ε, π , α) W s (ε, π , α, ∇α = 0), as,

R(α) = ε ∈ M s : ∂ W s ∂α ≥ 0 , R * (α) = σ ∈ M s : ∂ W * s ∂α ≤ 0 , ( . )
where M s represents the space of symmetric second-order tensors. W * s is the conjugate of W s obtained employing a Legendre-Fenchel transformation about the pair (σ , ε),

W * s (σ , π , α) = sup ε ∈M {σ . . ε -W s (ε, π , α)}. ( . )
Utilizing the expression of W s (ε, π , α) derived from Eq.( . ) the transformation Eq.( .) can be identi ed as,

W * s (σ , π , α) = 1 a(α) 1 2 σ . . S . . σ + btr(S . . σ )∆π + 1 2 b 2 tr(S . . I )∆π 2 + 1 2N ∆π 2 + ϕ 0 π -w(α), ( . )
and the elastic domains Eq.( .) can be obtained explicitly,

R(α) = ε ∈ M s : (ε -ε 0 ) . . C . . (ε -ε 0 ) ≤ w 1 1 -α , R * (α) = σ ∈ M s : σ . . S . . σ + 2b tr(S . . σ )∆π + b 2 tr(S . . I )∆π 2 ≤ w 1 (1 -α) 3 , ( . )
As α grows from 0, it is clear that the elastic domain in the total stress space reduces. is indicates a total stress-so ening behavior of the model as damage evolves. As well, it is interesting to note the dependency of the elastic domain in the total stress space on the equivalent pore pressure, π , and consequently on the saturation degree of the we ing uid, S w . is is because the total stress, σ (ε, π ), as mentioned earlier, involves a negative contribution due to the equivalent pore pressure, along with the e ective stress, σ e (ε), which is the same as that of an ordinary elastic solid. is dependency on the equivalent pore pressure can be recast into a dependency on the solution of the uid ow problem which is discussed further.

. . Transient hydraulic problem

e hydraulic part of the problem involves the spatio-temporal evolution of the uids within the porous skeleton under the in uence of gradients within uid pressures and imposed uxes at the boundaries. Bulk force due to gravity is neglected for this purpose. Since, we assume a passive and rare ed air phase, the hydraulic problem reduces to the solution of the Richards' equation for deformable porous skeleton. e corresponding developments were detailed in Sec.( . . ).

In particular the governing PDE, Eq.( . ), for the evolution of saturation degree of we ing uid, S w (x, t), is recapitulated below,

∂(ϕS w ) ∂t + κ η w ∇ • k w (S w ) (∇p c (S w )) = 0, ( . )
wri en in the "mixed water content" form and in the absence of gravity.

In the current chapter a speci c choice is done in order to reduce the non-linearity of the otherwise highly non-linear equation above. is concerns the algorithm utilized to resolve it numerically. e negative capillary pressure is replaced by the we ing uid pressure, p w (x, t) = -p c (S w ), which is considered also as an unknown. See Sec.( . . ) for similar developments leading to the "head" form of the Richards' equation. Consequently, we arrive at the following PDE with two unknowns, S w (x, t) and p w (x, t),

∂(ϕS w ) ∂t - κ η w ∇ • k w (S w ) (∇p w ) = 0. ( . )
As detailed further in Sec.( . ), an iterative approach is taken by rst considering S w (x, t) as a known quantity and solving the above equation for p w (x, t), in tandem with the mechanical problem for displacements of the solid skeleton. While doing this k w (S w ) is also a known function through Eq.( . ) and ϕ(p w , S w ) is given by the constitutive law Eq.( . ).

Appropriate boundary conditions accompany the resolution of Eq.( .), which could be either an imposed ux representing a natural boundary condition or an imposed pressure representing an essential boundary condition. Once p w (x, t) is resolved, S w (x, t) is obtained algebraically through the inverse retention relation Eq.( . ) recapitulated below,

S w (p w ) = p -1 c (p w ), ( . )
e iterations are then repeated until convergence.

Remark: e state variable π of the problem related to the skeleton, see Eq.( . ), is obtained from the solution pair, (p w , S w ), of the dissipative uid ow problem. us at every time t > 0, the solution of the hydraulic problem acts as an instantaneous input to the poromechanical damage problem through the parametrization of π . Consequently the bulk energy density in Eq.( . ) can be rewri en as,

W s (ε, α, ∇α; π ) = G s (ε, α; π ) + w(α) + 1 2 w 1 2 ∇α • ∇α . ( . )
All the developments done earlier still hold true under such parametrization.

. . Evolution problem

Assuming that the transient hydraulic problem is resolved to obtain π , we now focus on the evolution problem for the pair, (u, α). is is wri en using a variational approach similar to the developments in P M ( a,b). e principles governing the evolution are elaborated for the current problem which now concerns the porous skeleton.

e variational approach is based on a total energy of the porous skeleton associated to admissible pair of states ( , β) ∈ C ×D. ese functional spaces are de ned as,

C = { ∈ H 1 (Ω) n d : = 0 on ∂Ω D }, D = {β ∈ H 1 (Ω) : 0 ≤ β < 1 in Ω}, ( . )
where H 1 (Ω) is the Sobolev space of square integrable functions whose weak gradients are also square integrable. ∂Ω D ⊆ ∂Ω where displacement, u, is imposed.

is de nition of total energy at each time t is given as the integral over the whole domain of the bulk energy density minus the potential of external forces. In the absence of volume forces and surface forces acting respectively on Ω and its boundary ∂Ω, the potential of external forces vanishes and the total energy at a given time t reads,

E t ( , β) = ∫ Ω W t s ( , β, ∇β; π t (p t , S t )) dx, ( . )
Parametrization w.r.t π t of E t ( , β) is understood implicitly.

Remark:

We assume that all the elds are su ciently smooth in time allowing for the following developments. Also we consider only the states of damage where α < 1 and the saturation degree S w > 0. is is done so that in the following analysis the total energy remains nite.

e variational approach developed in P M ( a,b) for the evolution in Ω of (u t , α t ) ∈ C × D for all t ≥ 0 is then governed by the three principles of irreversibility of damage, stability and energy balance that read:

(a) Irreversibility of damage: t → α t must be non-decreasing. Consequently the admissible states accessible from α t can be de ned as,

D t = {β ∈ H 1 (Ω) : α t ≤ β < 1 in Ω}. ( . ) (b)
Stability: e state (u t , α t ) must be directionally stable in the sense that for all ( , β) ∈ C × D t , there exists h > 0, such that,

∀h ∈ [0, h], E t (u t + h( -u t ), α t + h(β -α t )) ≥ E t (u t , α t ). ( . )
(c) Energy Balance: During the evolution t → (u t , α t ) the following energy balance must hold,

E t (u t , α t ) = E 0 (u 0 , α 0 ) ( . )
where (u 0 , α 0 ) denotes the state of the skeleton at time t = 0.

. . First-order stability & necessary conditions e stability condition Eq.( .) essentially means that in the total energy landscape, within the immediate neighborhood of the state (u t , α t ) all the admissible states are of either higher or equal energy compared to that of (u t , α t ). Assuming this neighborhood, governed by h, to be small, one can expand the perturbed energy and rewrite Eq.( . ) as,

hE t (u t , α t )( -u t , β -α t ) + h 2 2 E t (u t , α t )( -u t , β -α t ) 2 + o(h 2 ) ≥ 0.
( . )

E t and E t represent respectively the rst and second directional derivatives of E t further referred to as FDD and SDD respectively for compactness.

e representation E t (u t , α t )(•) 2 is to be understood as a shorthand for the quadratic form whereas the associated symmetric bi-linear form is represented E t (u t , α t ) •, • i.e., the application of E t (u t , α t ) to the pair of Pore pressure, p w , and Saturation degree, S w , of the we ing uid at time t are further denoted in this chapter by p t and S t respectively in order not to overload the notation. directions •, • . e directional derivatives in Eq.( . ) have the following forms in the general direction ( , β),

E t (u t , α t )( , β) = ∫ Ω ∂W t s ∂ε . . ε( ) + ∂W t s ∂α β + ∂W t s ∂(∇α) • ∇ β dx, ( . ) E t (u t , α t )( , β) 2 = ∫ Ω a(α t )ε( ) . . C . . ε( ) + 2a (α t )(ε t -ε 0 ) . . C . . ε( ) β + 1 2 a (α t )(ε t -ε 0 ) . . C . . (ε t -ε 0 ) β 2 + w 1 2 ∇ β • ∇ β dx .
( . )

In Eq.( . ) the partial derivatives of bulk energy density are functions of state (u t , α t ) given by Eq.( . ). Dividing Eq.( .) by h and passing to the limit h → 0 gives,

E t (u t , α t )( -u t , β -α t ) ≥ 0, ( . )
which is the so-called rst-order stability condition and can be viewed as characterizing stationarity of the state (u t , α t ). In Eq.( . ) testing with β = α t and noting that C is a linear space can obtain the variational (weak) form of the classical equilibrium condition,

∫ Ω ∂W t s ∂ε . . ε( -u t )dx = 0, ∀ ( -u t ) ∈ C, ( . )
where the total stress tensor can obviously be identi ed as,

σ t = ∂W t s ∂ε = a(α)C . . (ε t -ε 0 ) -b [π t -π 0 ] I . ( . )
Using Eq.( . ) in Eq.( . ) gives the variational form of the non-local damage criterion,

∫ Ω ∂W t s ∂α (β -α t ) + ∂W t s ∂(∇α) • ∇(β -α t ) dx ≥ 0, ∀ β ∈ D t ( . )
Employing classical localization arguments of calculus of variations one can obtain from Eq.( . ) and Eq.( . ) the following local (strong) forms respectively with corresponding boundary conditions,

∇ • σ t = 0 in Ω, σ t • n = 0 on ∂Ω\∂Ω D , ( .
)

1 2 a (α t )(ε t -ε 0 ) . . C . . (ε t -ε 0 ) + w (α t ) -w 1 2 ∆α t ≥ 0 in Ω, ∂α t ∂n ≥ 0 on ∂Ω, ( . )
where n denotes the outwards unit normal vector associated to part of the boundary wherever invoked. Note the absence of surface and volume forces in the equilibrium equation according to earlier assumption. Since we have assumed the evolution is smooth in time, taking a time derivative of the Energy balance leads to,

0 = d dt E t (u t , α t ) = ∫ Ω ∂W t s ∂ε . . ε( u t ) + ∂W t s ∂α α t + ∂W t s ∂(∇α) • ∇ α t dx . ( . )
Integrating by parts the gradient terms and further using the local form of the equilibrium equations Eq.( . ) gives,

0 = ∫ Ω ∂W t s ∂α -∇ • ∂W t s ∂(∇α) α t dx + ∫ ∂Ω ∂W t s ∂(∇α) • n α t dx . ( . )
Owing to the irreversibility of damage everywhere in Ω and the local inequalities Eq.( .), the two integrals on the right-hand side of the above equation are non-negative. So, both of them should vanish. Further using classical localization arguments we obtain the so-called consistency conditions or the Karush-Kuhn-Tucker (KKT) conditions applicable everywhere in Ω and on the boundary ∂Ω respectively,

1 2 a (α t )(ε t -ε 0 ) . . C . . (ε t -ε 0 ) + w (α t ) -w 1 2 ∆α t α t = 0 in Ω.
∂α t ∂n α t = 0 on ∂Ω.

( . )

ese conditions can be read in tandem with the Irreversibility of damage. e rst condition states that everywhere in Ω, damage increases only if the local form of the damage criterion Eq.( .) is an equality and if it is a strict inequality then damage does not increase. e second condition states that everywhere on the boundary, ∂Ω, if damage increases then the spatial derivative normal to the boundary vanishes.

. . Tension-Compression split

e regularized model for fracture with unilateral contact at the fracture faces proposed in A . ( ) is adopted here with an intention to avoid possible interpretation of the fracture faces. is is based on the decomposition of the strain tensor into its spherical and deviatoric components,

ε = 1 n d ϵI + ε D . ( . )
Further a decomposition of ϵ into positive, ϵ + = 0.5(ϵ + |ϵ |), and negative, ϵ -= 0.5(ϵ -|ϵ |), parts is done to distinguish the contributions within elastic energy density due to dilation and compression. A ecting the degradation due to damage only to the contributions due to dilation and shear is shown (A ., ) to avoid crack interpenetration by predicting asymmetric results under traction and compression tests. Such a split of energy density is adopted such that,

G s (ε, π , α) = a(α)G + s (ε) + G - s (ε, π ), G + s (ε) = 1 2 Kϵ 2 + + µε D . . ε D , G - s (ε, π ) = 1 2 Kϵ 2 --bϵ∆π - 1 2N ∆π 2 -ϕ 0 π , ( . )
in the absence of pre-strain. K and µ are respectively the bulk and shear moduli of the intact skeleton.

While the earlier developments still remain valid, the above split amounts to two modi cations. e rst one concerns the state equations. Speci cally, Eqs.( . (a) & (c)) are modi ed to,

∂W s ∂ε = a(α) (Kϵ + I + 2µε D ) + Kϵ -I -b∆πI , ∂W s ∂α = a (α) 1 2 Kϵ 2 + + µε D . . ε D + w (α). ( . )
e total stress tensor and the e ective stress tensor now read respectively as,

σ = σ e -b∆πI , σ e = a(α) (Kϵ + I + 2µε D ) + Kϵ -I . ( . )
e resulting elastic domains within the stress and strain spaces, Eq.( .), are modi ed to,

R(α) = ε ∈ M s : Kϵ 2 + + 2µε D . . ε D ≤ w 1 1 -α , R * (α) = tr + (σ ) 2 n 2 d K + 1 2µ σ D . . σ D + 2b n d K [tr -(σ )∆π + + tr + (σ )∆π ] + b 2 K ∆π 2 + ≤ w 1 (1 -α) 3 , ( . )
where ∆π ± = 0.5(∆π ± |∆π |) and tr ± (σ ) = 0.5(tr(σ ) ± |tr(σ )|). e second modi cation concerns the damage criterion, Eq.( . ), and the KKT conditions in the bulk, Eq.( . (a)). Due to the above split these conditions are now wri en respectively as,

a (α) 1 2 Kϵ 2 + + µε D . . ε D + w (α t ) -w 1 2 ∆α t ≥ 0 in Ω, , ( . )
a (α) 1 2 Kϵ 2 + + µε D . . ε D + w (α t ) -w 1 2 ∆α t α t = 0 in Ω. ( . )
. .

Fracture toughness & related assumptions

Fracture in the context of gradient damage models is understood as the localization of damage over a nite length controlled by the internal length, . In the works of P M ( ); M . ( ) a link between fracture toughness, G c , of a material and damage localization has been derived by explicitly constructing the total energy of a non-homogeneous damage evolution in one-dimension and passing to the limit of zero stress within a localized zone representative of a fracture. e energy dissipated in the creation of a single fracture is then identi ed with G c . e following relation is thus arrived at in a more general se ing of damage laws,

G c = c w w 1 √ 2 with c w = 4 1 ∫ 0 w(β) w 1 dβ . ( . )
For the particular forms of the damage laws chosen, Eqs.( . & . ), the above relation reads,

G c = 4 √ 2 3 w 1 , c w = 8 3 . ( . )
A similar analysis in the current context would lead to an identi cation of G c which is not just dependent on the damage laws chosen but also on the saturation degree through the dependency on equivalent pore pressure. Such a dependency was proposed in the work of C . ( ). In fact, in the experimental work of L . ( ) it has been observed that fracture toughness of soils decreases exponentially with increasing water content. Also, in a uni-axial tension test one can infer from Eq.( . ) that the elastic limit stress, σ e , beyond which damage starts to initiate is linked to π .

Algorithm : Alternating algorithm for poromechanical gradient damage model

Input:

(u n-1 , π n-1 , α n-1 ) initialization: (u 0 n , π 0 n , α 0 n ) ← (u n-1 , π n-1 , α n-1 ) ; i = 0 ;
while alternate algorithm not converged do • initialization:

(u i n 0 , π i n 0 ) ← (u i-1 n , π i-1 n ) ; • k = 0 ; while poro-split algorithm not converged | α n = α i-1 n do * Solve (u i n k , p i n k ) {Equlibrium, Eq.( . ) & pressure evolution, Eq.( . )}; * Compute S i n k (p i n k ) {Inverse retention relation, Eq.( . )}; * Compute π i n k (p i n k , S i n k ) {Equivalent pore pressure, Eq.( . )}; * k ← k + 1; end • update: (u i n , π i n ) ← (u i n k , π i n k ); • Solve α i n argmin α E n (u i n , α) | α n-1 ≤ α ≤ 1 in Ω; • i ← i + 1; end update: (u n , π n , α n ) ← (u i n , π i n , α i n ); Output: (u n , π n , α n )
e minimization problem for damage within each alternate iteration is posed as,

α i n = argmin α E n (u i n , α) | α n-1 ≤ α ≤ 1 in Ω, ( . )
where the unilateral constraint α n-1 ≤ α ≤ 1 is the time-discrete version of the irreversibility of damage. e spatial discretization is done employing linear Lagrange elements with characteristic element size equal to half of the internal length, ∆x = 0.5 . e above minimization problem is resolved using a bound-constrained optimization solver routine available as part of the TAO library (B .,

). e convergence criterion of the alternating algorithm, at each iteration i, is the comparison against a tolerance, the 2norm of the di erence between the damage solutions of successive iterations,

(α i n -α i-1 n ) .
.

Table . : Material properties, model parameters used through Chp.( ), unless mentioned otherwise.

property/ parameter

Clay

ϕ [-] . E [Pa] . E ν [-] . b [-] . N [Pa] . E κ [m 2 ] . E- η w [Pa.s] . E- π 0 [Pa] m [-] . S r es r [-] . G c [N.m -1 ] . [m] . k [-]
. E-e above modeling approach is now applied to study the desiccation of soils. Typical desiccation experiments are done by subjecting a certain mass of fully saturated soil to air drying under controlled temperature and relative humidity. A few test cases of desiccation experiments can be found in the literature, P . ( ); S S ( , a); S ( ) to name a few. In what follows, we intend to verify if the modeling approach above results in localization of damage solutions thus allowing fracturing in drying soils. In order to de nitely verify this, one would need to consider homogeneous solutions, at least along one direction and look for bifurcations leading to localizations in that particular direction.

A test case similar in construction is that of a thermal shock as studied by S . ( ) in which homogeneous damage along the surface of a ceramic plunged into a cold bath has been studied for possible bifurcations as time passes. We take inspiration from this study, to impose boundary conditions, as detailed in Sec.( . . ), such that drying of an initially saturated sample occurs homogeneously along the drying surface at initial times.

Further we consider a plane-strain assumption owing to the transverse dimension of the samples in all the experiments being larger than the vertical depth. is assumption does not a ect drastically the above developments. Speci cally, the in-plane stress components can still be obtained from the relation Eq.( . ) with the de nition of sti ness tensor in index notation,

C i jkl = λδ i j δ kl + µ(δ ik δ jl + δ il δ jk ), i, j, k, l ∈ {1, 2}.
( . )

λ and µ are respectively the rst and second Lamé parameters related to the empty porous skeleton. So, the boundary value problem formed by the coupled system of equations Eqs.( . , . ) and the bound-constrained minimization w.r.t α, Eq.( .), are resolved with appropriately de ned boundary condition and initial conditions as laid out further in Sec.( . . ) and following the algorithm described in Alg.( ). e material properties of the porous medium and the parameters of the model chosen for the purpose of the simulations are listed in Table .( . ), which are in the range typical of clays saturated with air-water mixture.

. .

Problem setup

= {(x 1 = 0)∪(x 1 = L)∪(x 2 = 0) ∪ (x 2 = H )}. At time t = 0 it
is assumed that the porous skeleton is completely intact, stress free and fully saturated giving, α 0 (x) = 0, u 0 (x) = 0, ε 0 (x) = 0, π 0 (x) = 0.

( . )

x ∈ Ω is the position vector given by x = x 1 e 1 + x 2 e 2 . e horizontal and vertical displacements are set to vanish respectively on the lateral (x 1 = 0, x 1 = L) and the bo om (x 2 = H ) faces along with the shear stresses at those boundaries. e top face at x 2 = 0 is a free surface. ese set of mechanical boundary conditions ∀t ≥ 0 read,

u t • e 1 (x 1 =0) = 0, u t • e 1 (x 1 =L) = 0, u t • e 2 (x 2 =H ) = 0, e 2 • σ t • (-e 1 ) (x 1 =0) = 0, e 2 • σ t • e 1 (x 1 =L) = 0, σ t • (-e 2 ) (x 2 =0) = 0, e 1 • σ t • e 2 (x 2 =H ) = 0.
( . )

Damage is not prescribed and is free to evolve at all the boundaries whenever the damage criterion is met. Accordingly the natural boundary condition on damage reads,

∂α t ∂n x ∈∂Ω = 0 ∀t ≥ 0. ( . )
For the hydraulic problem, the lateral and bo om faces are set to be impermeable. e loading is done by a constant imposed suction, p < 0, on the top face thus inducing a drying e ect. ese boundary conditions ∀t ≥ 0 read,

∇p t • (-e 1 ) (x 1 =0) = 0, ∇p t • e 1 (x 1 =L) = 0, p (x 2 =0) = p, ∇p t • e 2 (x 2 =H ) = 0.
( . )

. . Results

It can be seen from the problem setup that the geometry, material properties and initial conditions are invariant along x 1 -direction. e loading and boundary conditions on the top and bo om boundaries are as well invariant along x 1 -direction. e boundary conditions on the lateral faces demand vanishing gradients of the solution along the x 1 -direction. So, even though the problem does not render itself for obtaining an easy exact solution due to its non-linearity, one can a priori anticipate the existence of a particular class of solutions that are homogeneous in the x 1 -direction and are dependent only on x 2 . Bifurcation from such solutions would mean that at some time t > 0 a dependency of the solution on x 1 arises. is leads to localization within these pro les with a so-called period doubling phenomenon observed in similar applications of gradient damage modeling (B ., ; S ., ; M ., ). In essence selective initiation and propagation of fractures occur resulting a complex fracture network in higher dimensions. ese localized fractures propagate further into the domain as the desaturation progresses.

e periodicity and spacing is qualitatively similar to what was observed in experimental desiccation tests as noted in Sec.( . . ).

When it comes to the driving force behind the damage growth in this problem setup, one needs to look at the x 1 -homogeneous solution structure prior to any bifurcation. Fig. ( . ) shows the solution pro les along the depth of the domain at time t = 0.8s. While the pressure and the saturation pro les at the top boundary, x 2 = 0m, con rm to the boundary condition imposed, the solution within the pro le is only a ected up to a certain depth controlled by the uid ow problem and this depth grows in time. e compression induced due to this desaturation is apparent looking at the pro le of ϵ -, whereas ϵ + remains vanishing. ere is as well a signi cant contribution due to the deviatoric strain, ε D , which in fact helps to drive the damage in accordance with Eq.( . ). is can be seen in e ow of we ing uid within the domain is governed by the gradient of pressure and the saturation degree according to the Darcy's law, Eq.( . ). Fig. ( . ), shows both the components of the volumetric uid ow vector, V evolution along x 1 at various depths when the series of fractures have already localized (approximate locations plo ed as red vertical lines). It is clear that the vertical component of ow, V • e 2 , is not signi cantly a ected due to the presence of a fracture within the domain. is is obviously not in line with the prevailing knowledge that uid ows more freely along the fracture plane due to its higher permeability compared to intact porous solid. is is a direct consequence of the absence of coupling in the current modeling approach, between the damage eld and the permeability of the material. e horizontal component of the ow, V • e 1 , indicates a tendency of asymmetric ow in the vicinity of the fracture tip where uid seems to ow into the tip ahead of it and out of the tip in its At the top boundary and along the length of the fracture there is not a signi cant e ect.

. e presented modeling approach is similar in construction to the prevailing paradigm (M M , ; C ., ; H S ,

) that the creation of opening mode fracture in desiccating soils is due to the release of elastic energy. Apart from the two-dimensional test case presented, one can imagine as well partially constrained (bo om boundary xed, lateral boundaries free) desiccation test cases which could generate also tensile e ective stresses on the drying face as observed in the aforementioned literature. e choice of boundary conditions in the current study allowed us to con rm that x 1 -homogeneous solutions could bifurcate in a periodic manner and localize as fractures. e intent of this chapter is to investigate a phenomenon of fracture initiation driven by capillarity, especially applicable to ne-grained soils close to full saturation. In order to do this in the framework of macroscopic porous continuum, the required ingredient is postulated as a damage mechanism described by a scalar damage variable, much like the gradient damage models that describe fractures in bri le materials. However, the driving force behind the damage evolution is revisited to replicate the capillary forces that are in action at the pore-scale. e analysis starts with the postulation of a damaged porous solid and its associated energy density. is postulate is a phenomenological one inspired by experimental observations and is designed to be simplistic in order to allow a detailed analysis. Following this, the evolution problem for damage is derived starting from a total energy of the porous solid and in accordance with the variational approach based on the three principles of irreversibility of damage, stability and energy balance, assuming a minimization principle to hold true at each time. is lays the foundation for investigation of the possibility of fracture initiation in a two-dimensional desiccation problem.

Similar to the previous chapter, the phase eld model of partial saturation is not adopted.

e starting points of the current development are the classical poromechanical approach to partial saturation (see Sec.( . )), and the variational approach underlying the gradient damage modeling of bri le fracture (see M . ( ) for an overview). In addition to this we work under an overarching assumption of a hydrophilic porous skeleton whose pore spaces are lled by two uids, air and water, between which air is considered passive and the hypothesis of small perturbations is assumed as detailed in Sec.( . . ). Damage of the porous solid is assumed to induce a reduction in its airentry pressure and in particular to degrade its retention properties. Since according to the classical unsaturated poromechanics (see C ( )), a porous solid is understood as a so-called 'we ed' porous skeleton with a thin layer of uid a ached to the pore walls, the aforementioned degradation is equivalent to reduction of the cumulative interfacial energy stored at the interfaces between the uids and the solid pore walls as well as at the uid-uid interfaces.

Accordingly the degraded retention relation between the saturation degree of the we ing uid, S w and the degraded macroscopic capillary pressure, p c , reads,

p c (S w , α) = a(α) p c (S w ) = -a(α) p w (S w ) = -p w (S w , α), ( .
)
where α plays the role of a scalar damage variable and a(α) is a damage law whose functional form is elaborated at a later point. With an abuse of notation p c (S w ) now represents the standard retention relation of the form Eq.( . ). Also, the assumptions leading to Richard's equation are adopted, in particular the hypothesis of a passive non-we ing phase is assumed, leading to p c ≈ -p w , with p w as the pore-water pressure in a degraded porous solid and p w as its non-degraded counterpart. While the above serves as a reasonable rst assumption for investigation, a complete modeling approach should also involve the degradation of the resistance to uid ow within damaged zone and eventually along the localized fracture plane. We consider again an n d -dimensional domain, Ω ∈ n d , representing the initial con guration of a damaging isotropic partially saturated porous solid whose boundary is denoted by ∂Ω. In line with standard gradient damage modeling (M ., ) and the theory of partially saturated poroelasticity, the following considerations are done:

. e scalar internal variable, α ∈ [0, 1], is assigned the role of describing the extent of damage in the sense of loss of retention properties. α = 0 denotes an intact healthy skeleton whose retention properties are described by the standard retention curve. Whereas, α = 1 denotes a fully damaged skeleton whose retention properties are degraded and can only hold vanishing or residual amount of we ing uid. . For a given unit volume of porous solid in its reference con guration the bulk energy density, W , accounting for the internal energy and the dissipation due to material degradation, is a state function characterized by: (ε, ϕ, S w , α, ∇α), respectively the linearized strain tensor, the Lagrangian porosity, the saturation degree of the we ing uid, the damage variable and the gradient of damage.

. In a similar way as the construction of the standard gradient damage model (M ., ) for the solid continuum, the porous solid is assumed to be dissipative in a non-local sense due to the dependency of the bulk energy density on the gradient of damage. e particular expression of the bulk energy density is assumed to be:

W (ε, ϕ, S w , α, ∇α) = Ψ s (ε, ϕ, S w , α) + w(α) + 1 2 w 1 2 ∇α • ∇α . ( . )
e di erent terms of the above expression are explained below.

(a) e energy density of the porous solid, Ψ s (ε, ϕ, S w , α), as elaborated in Ch.( ), encompasses the contributions due to the deformation of the porous skeleton, resulting in solid strains and changes in porosity and also, the interfacial energy contribution due to formation/annihilation of interfaces. Adopting the concept of energy separation (C , ) between the free energy of the porous skeleton and the interfacial energy we obtain, Ψ s (ε, ϕ, S w , α) = ψ s (ε, ϕ) + ϕ U (S w , α).

( . )

Here a speci c choice is done that would result in the aforementioned postulate that damage means the degradation of only the retention properties and not the poroelastic properties. e la er was the choice for modeling tensile mode fracture in Ch.( ). e choice done here allows to isolate the investigation of fracture formation just due to the 'release' of interfacial energy. Specically, we assume a degradation of the form,

U (S w , α) = a(α) U (S w ), ( . )
consistent with the de nition of the degraded capillary pressure in Eq.( . ). Again with an abuse of notation U (S w ) now represents the undamaged interfacial energy given by,

U (S w ) = 1 ∫ S w p c (s)ds. ( . )
While various choices of the functional form of a(α) are possible, we make a speci c choice following the standard gradient damage modeling,

a(α) = (1 -α) 2 + k , ( .
)
where k is a small positive constant used solely for numerical purposes to govern the fully damaged state. Accordingly, when α grows from 0 to 1, the interfacial energy density degrades to a vanishing value. e behavior of the capillary pressure and interfacial energy density can be visualized in the Fig. ( . ). e consequence on the damage evolution due to the above choices is discussed further. Now concerning the free energy density of the skeleton, a possible form that retrieves back the classical constitutive relations as its state equations reads as,

ψ s (ε, ϕ) = 1 2 (ε -ε 0 ) . . C . . (ε -ε 0 ) + 1 2 b 2 N (ϵ -ϵ 0 ) 2 -∆ϕ(bN (ϵ -ϵ 0 ) -p * 0 ) + 1 2 N ∆ϕ 2 , ( .
)
where the notations are borrowed from Ch.( ). ∆ϕ represents the increment of Lagrangian porosity, (ϕ -ϕ 0 ). ϕ 0 , ε 0 and p * 0 denote the initial reference states of the Lagrangian porosity, the linearized strain tensor and the average pore-pressure respectively.

(b) e second term in Eq.( .) is the local part of the dissipated energy,

w(α) = w 1 α, ( . )
growing from 0 when α = 0, to a positive constant w 1 < +∞ when α = 1. us in accordance with the developments in P . ( a), since w (α) > 0 there exists an elastic phase preceding damage initiation and the niteness of w 1 ensures that the energy dissipated during a homogeneous evolution of α from 0 to 1 is as well nite. However, w 1 is not related to the fracture toughness in the classical sense. Rather it is to be viewed as related to a threshold beyond which the grains that form the skeleton start to 'slip' at their contacts thus characterizing a dissipative phenomenon that accompanies the release of built-up interfacial energy.

(c) e last term in Eq.( .) is the non-local dissipation, which is assumed to be a quadratic function of the gradient of α, and is intended to regularize the local model allowing localizations of nite thickness. appears with the physical dimension of a length that in the context of gradient damage modeling is intended to control the localization thickness.

. e dual relations associated to the state variables in Eq.( .) are obtained as follows,

∂W ∂ε = C . . (ε -ε 0 ) + bN (b(ϵ -ϵ 0 ) -∆ϕ)I , ∂W ∂ϕ = -bN (ϵ -ϵ 0 ) + p * 0 + N ∆ϕ + U (S w , α), ∂W ∂S w = ϕ ∂U ∂S w (S w , α), ∂W ∂α = ϕa (α) U (S w ) + w (α),

∂W ∂(∇α)

= w 1 2 ∇α .

( . )

. .

Transient hydraulic problem

Similar to the what was said in Sec.( . . ), the transient hydraulic problem governing the evolution the uids within the porous solid is considered with the absence of bulk force due to gravity. is evolution can be determined by resolving Richards' equation for deformable porous skeleton, Eq.( . ), that governs the pore pressure of the we ing uid, p w , recapitulated below for easy reference,

∂(ϕS w ) ∂p w ∂p w ∂t + κ η w ∇ • kw (p w ) (-∇p w ) = 0. ( . )
Unlike the developments in Sec.( . . ), S w is now obtained from the degraded inverse retention relation,

S w (p w , α) = p -1 c -p w a(α) , ( .
)
where p -1 c (•) is the undamaged inverse retention relation. e intrinsic permeability κ is assumed independent of α. However, since κ is expected to increase, in orders of magnitude, within a fracture we make a rather simplifying assumption of kw (p w ) = 1 to in part compensate for the lack of the aforementioned coupling. While this serves as a reasonable rst assumption for investigation, a complete modeling approach should also involve the degradation of the resistance to uid ow within damaged zone and eventually along localized fracture planes.

Appropriate boundary conditions compliment the above PDE which could be either an imposed ux representing a natural boundary condition or an imposed pressure representing a Dirichlet boundary condition.

Remark: u, the displacement of the solid, ϕ and α which are the solutions of the poromechanical damage problem are thus coupled to the hydraulic problem at each instant. Assuming a possible split between the resolution of both problems at each instant, the solution of the hydraulic problem can be viewed as providing instantaneous input to the poromechanical damage problem through the parametrization of S w . Consequently the bulk energy density in Eq.( .) can be rewri en as,

W (ε, ϕ, α, ∇α; S w ) = Ψ s (ε, ϕ, α; S w ) + w(α) + 1 2 w 1 2 ∇α • ∇α . ( . )
All the developments done earlier still hold true under such parametrization.

. . Evolution problem

Now we pose again the evolution problem for the solution triplet (u, ϕ, α) of the poromechanical damage problem using the variational approach similar to the developments in Sec.( . . ). At the risk of reintroducing this approach, the principles governing the evolution are elaborated for the current problem which now concerns the porous solid and not just the skeleton. e variational approach starts again with the de nition of a suitable total energy of the body under consideration (i.e., the porous solid), that is associated to the triplet of admissible states ( , φ, β) ∈ C × P × D. ese functional spaces C and D are as de ned in Eq.( . ) and P is de ned as,

P = {φ ∈ L 2 (Ω) : 0 < φ < 1 in Ω}, ( .
)
where L 2 (Ω) is the Lebesgue space of square integrable functions equipped with an L 2 -norm. ∂Ω D ⊆ ∂Ω where displacement, u, is imposed. is de nition of total energy at each time t is given as the integral over the whole domain of the bulk energy density minus the work of external forces,

E t ( , φ, β) = ∫ Ω W t ( , φ, β, ∇β; S t ) dx -W e t . ( . )
Since the focus here is on the porous solid, the external loading concerns not only bulk forces and tractions exerted on the solid skeleton, but also zeroth-order bulk actions tending to modify from inside the porosity of the solid. is time dependent bulk action is obviously related to the presence of uids in the porous network and is therefore coupled with the solution of the hydraulic problem. In a similar way as in the case of fully saturated porous materials, working this action of the change of Lagrangian porosity it will coincide with the average pore pressure , p * t = p t S t + p nw S nw , for all t > 0. e assumption of a passive non-we ing phase leads to p * t (S t , α t ) ≈ p t S t . Note the dependence of p * t on damage owing to the coupling of the uid problem to the damage evolution within the domain. is contribution to the potential due to external e orts de ned on the set of admissible porosity Pore pressure, p w , and Saturation degree, S w , of the we ing uid at time t are further denoted in this chapter by p t and S t respectively in order not to overload the notation.

elds, φ ∈ P, reads,

W e t (φ) = ∫ Ω p * t φdx . ( . )
Remark: We assume that all the elds are su ciently smooth in time allowing for the following developments. Also we consider only the states of damage where α < 1 and the saturation degree S w > 0. is is done so that in the following analysis the total energy remains nite.

Now we are in a position to setup the principles of variational approach for the evolution in Ω of (u t , ϕ t , α t ) ∈ C × P × D for all t ≥ 0 which read:

(a) Irreversibility of damage: t → α t must be non-decreasing. Consequently the admissible states accessible from α t belong to the space D t de ned in Eq.( . ) (b) Stability: e state (u t , ϕ t , α t ) must be directionally stable in the sense that for all ( , φ, β) ∈ C × P × D t , there exists h > 0, such that,

∀h ∈ [0, h], E t (u t + h( -u t ), ϕ t + h(φ -ϕ t ), α t + h(β -α t )) ≥ E t (u t , ϕ t , α t ). ( . )
(c) Energy balance: During the evolution t → (u t , ϕ t , α t ) the following energy balance must hold,

E t (u t , ϕ t , α t ) = E 0 (u 0 , ϕ 0 , α 0 ) - t ∫ 0 ∫ Ω ϕ s ∂π ∂S s S s dx ds, ( . )
where the classical de nition of equivalent pore-pressure has been adapted to the current context, π (S t , α t ) = p * t (S t , α t ) -U (S t , α t ), and (u 0 , ϕ 0 , α 0 ) denotes the state of the skeleton at time t = 0. At any time 0 ≤ s ≤ t, the second term on the le hand side of Eq. . accounts for the time parametrization of the total energy through its dependency on S t .

. . First-order stability & necessary conditions

As seen in Sec.( . . ), equilibrium equation and damage criterion can be obtained as rst-order stability conditions starting from Eq.( . ). In addition, since the solution now involves ϕ t , an associated zeroth-order balance equation is expected to appear. We start again by expanding the perturbed energy up to the second-order,

hE t (u t , ϕ t , α t )( -u t , φ -ϕ t , β -α t ) + h 2 2 E t (u t , ϕ t , α t )( -u t , φ -ϕ t , β -α t ) 2 + o(h 2 ) ≥ 0.
( . )

E t and E t represent respectively the rst and second directional derivatives of E t further referred to as FDD and SDD respectively for compactness. e representation E t (u t , ϕ t , α t )(•) 2 is to be understood as a shorthand for the quadratic form whereas the associated symmetric bi-linear form is represented E t (u t , ϕ t , α t ) •, • i.e., the application of E t (u t , ϕ t , α t ) to the pair of directions •, • . e directional derivatives in Eq.( . ) have the following forms in the general direction ( , φ, β),

E t (u t , ϕ t , α t )( , φ, β) = ∫ Ω ∂W t ∂ε . . ε( ) + ∂W t ∂ϕ -p * t φ + ∂W t ∂α -p * t ϕ t β + ∂W t ∂(∇α) • ∇ β dx, ( . ) E t (u t , ϕ t , α t )( , φ, β) 2 = ∫ Ω ε( ) . . C . . ε( ) + N (bϵ( ) -φ) 2 -ϕ t π t β 2 -2π t φ β + w 1 2 ∇ β • ∇ β dx, ( . )
where for the sake of compactness of the notation we denoted the directions ( -u t , φϕ t , βα t ) with ( , φ, β) and de ne their associated function spaces as C × P × D with,

P = { φ ∈ H 1 (Ω) : -1 < φ < 1 in Ω}. ( . )
In Eq.( . ) the partial derivatives of bulk energy density are functions of state (u t , ϕ t , α t ) given by Eq.( . ). Dividing Eq.( .) by h and passing to the limit h → 0 gives,

E t (u t , ϕ t , α t )( -u t , φ -ϕ t , β -α t ) ≥ 0, ( . )
which is the so-called rst-order stability condition and can be viewed as characterizing stationarity of the state (u t , ϕ t , α t ). In Eq.( . ) testing with φ = ϕ t , β = α t and noting that C is a linear space one can obtain the variational (weak) form of the classical equilibrium condition,

∫ Ω ∂W t ∂ε . . ε( -u t )dx = 0, ∀ ( -u t ) ∈ C, ( . )
where a stress tensor can obviously be identi ed as,

σ t = ∂W t ∂ε = C . . (ε t -ε 0 ) + bN (b(ϵ t -ϵ 0 ) -∆ϕ t )I . ( . )
Similarly testing with = u t , β = α t in Eq.( . ) one obtains a zeroth-order balance equation associated to small variations in ϕ t ,

∫ Ω ∂W t ∂ϕ -p * t (φ -ϕ t )dx = 0, ∀ φ ∈ P, ( . )
Finally, using Eq.( . ) and Eq.( . ) in Eq.( . ) gives the variational form of the non-local damage criterion,

∫ Ω ∂W t ∂α -p * t ϕ t (β -α t ) + ∂W t ∂(∇α) • ∇(β -α t ) dx ≥ 0, ∀ β ∈ D t ( . )
Employing classical localization arguments of calculus of variations one can obtain from Eq.( . ), Eq.( . ) and Eq.( . ) the following local (strong) forms respectively with corresponding boundary conditions, ) where n denotes the outwards unit normal vector associated to part of the boundary wherever invoked. Note the absence of surface and volume forces in the equilibrium equation according to earlier assumption. Eq.( .) is the local form of the zeroth-order balance law associated to variations in porosity. Rearranging it one can identify the relation which in classical poromechanics, Eq.( . ), is characterized as the constitutive relation for porosity,

∇ • σ t = 0 in Ω, σ t • n = 0 on ∂Ω\∂Ω D , ( . ) -bN (ϵ t -ϵ 0 ) + p * 0 + N ∆ϕ t + U (α t ; S t ) = p * t in Ω, ( . ) -π t ϕ t + w (α t ) -w 1 2 ∆α t ≥ 0 in Ω, ∂α t ∂n ≥ 0 on ∂Ω, ( .
ϕ t = b(ϵ t -ϵ 0 ) + 1 N [(p * t -p * 0 ) -U (α t ; S t )] + ϕ 0 . ( . )
Remark: e above formulation results in an equivalent equilibrium equation to be resolved in tandem with the zeroth-order balance law for variations in porosity, compared to the formulation in Chp.( ). is can be seen by substituting Eq.( . ) into the equilibrium equation, Eq.( . ), giving,

∇ • C . . (ε -ε 0 ) -b[(p * t -p * 0 ) -U (α t ; S t )]I = 0 in Ω, ( . )
which is similar in structure compared to Eq.( . ). In the above the classical stress tensor for unsaturated poroelasticity can be identi ed as, ) is equivalence between the formulations is exploited further for the purpose of numerical approximation.

σ t = C . . (ε t -ε 0 ) -b[(p * t -p * 0 ) -U (α t ; S t )]I . ( .
Since we have assumed the evolution is smooth in time, taking a time derivative of the Energy balance leads to,

0 = d dt E t (u t , ϕ t , α t ) + ∫ Ω ϕ t ∂π t ∂S t S t dx = ∫ Ω ∂W t ∂ε . . ε( u t ) + ∂W t ∂ϕ -p * t ϕ t + ∂W t ∂α -p * t ϕ t α t + ∂W t ∂(∇α) • ∇ α t dx . ( . )
Integrating by parts the gradient terms and further using the local form of the equilibrium equations Eq.( . ) and the zeroth-order balance law Eq.( . ) gives,

0 = ∫ Ω ∂W t ∂α -p * t ϕ t -∇ • ∂W t ∂(∇α) α t dx + ∫ ∂Ω ∂W t ∂(∇α) • n α t dx . ( . )
Owing to the irreversibility of damage everywhere in Ω and the local inequalities Eq.( .), the two integrals on the right-hand side of the above equation are non-negative. So, both of them should vanish. Further using classical localization arguments we obtain the KKT conditions applicable everywhere in Ω and on the boundary ∂Ω respectively,

-π t ϕ t + w (α t ) -w 1 2 ∆α t α t = 0 in Ω.
∂α t ∂n α t = 0 on ∂Ω.

( . )

ese conditions can be read in tandem with the Irreversibility of damage. e rst condition states that everywhere in Ω, damage increases only if the local form of the damage criterion Eq.( .) is an equality and if it is a strict inequality then damage does not increase. e second condition states that everywhere on the boundary, ∂Ω, if damage increases then the spatial derivative normal to the boundary vanishes.

.

&

In Chp.( ) an alternating algorithm inspired from the alternate minimization (B ., ) used extensively for the standard gradient damage modeling is described, Alg.( ). As noted in the previous chapter, the alternate minimization algorithm is based on the observation that while the minimization problem of the regularized energy functional is non-convex, when minimized separately for u and α, the individual minimization problems are convex. e same can be said to be true in the current context for minimizations w.r.t u, ϕ and α.

If one assumes the absence of damage modeling, numerical approximation of the poromechanical problem can be done using the nite element method and resolving in tandem the equilibrium equation, Eq.( .) and the zeroth-order balance law for variations in porosity, Eq.( . ), in their variational forms subject to prescribed boundary conditions. As per earlier remark, here we exploit the equivalence between the aforementioned formulation and the classical equilibrium equation of unsaturated poroelasticity, Eq.( . ). Speci cally we resolve the la er for the poromechanical part of the problem. For the transient hydraulic problem on the other hand we resolve the mixed-head form of the governing equation for pore-water pressure, Eq.( . ), with Eq.( . ). In the presence of damage, which is the current context, the damage problem is resolved by the minimization of the total energy under the unilateral constraint of irreversibility of damage.

Given the solution triplet (u n-1 , p n-1 , α n-1 ) of the solid and uid problems at time-step (n -1), Alg.( ) describes the alternating algorithm to obtain the solution at time-step n. Spatial discretization is done employing the Taylor-Hood elements for approximating the pair (u, p) and linear Lagrange nite elements for α. e time derivative within Eq.( . ) for the hydraulic problem is discretized using the implicit Euler scheme of rst-order.

e coupled problem for solid displacement and uid pressure, (u i n , p i n ), at each alternate iteration i is solved employing Newton iterations and a sparse LU decomposition routine, available in the FEniCS suite (A ., ), to solve the linearized systems. e minimization problem for damage within each alternate iteration is posed as, ) where the unilateral constraint α n-1 ≤ α ≤ 1 is the time-discrete version of the irreversibility of damage. Numerically we solve this minimization problem using a bound-constrained optimization solver routine available as part of the TAO library (B .,

α i n = argmin α E n (u i n , ϕ i n , α) | α n-1 ≤ α ≤ 1 in Ω, ( .
). e convergence criterion of the alternating algorithm, at each iteration i, is the comparison against a tolerance, the 2 -norm of the di erence between the damage solutions of successive iterations, (α i n -α i-1 n ) .

[-] . E [Pa] . E ν [-] . b [-] . N [Pa] . E κ [m 2 ] . E- η w [Pa.s] . E- π 0 [Pa] . m [-] . S r es r [-] . w 1 [N.m -3 ] [m] .
During the experiments in the aforementioned literature, two stages were observed in the drying process. e rst stage involving large irreversible deformations with the degree of saturation close to and the second stage involving noticeable decrease in the saturation degree and with smaller deformations. Fracture initiation usually has been associated to the sample close to full saturation and the air-water interface coinciding with the apparent soil surface. In the current modeling context damage initiation is associated to the threshold, w 1 , appearing within the local dissipation contribution, Eq.( . ), to the bulk energy density and this, as mentioned earlier, is understood as not directly related to a fracture toughness of the material but to the creation of a new uid-uid interface within the porous medium. Nevertheless, w 1 is a parameter that depends on the material and boundary conditions. For instance, H . ( ) have shown experimentally that the fracturing is the preferred mode of invasion of air into water saturated granular media when the con ning stress is lower. So, this threshold can only be characterized through experimental observation of fracture initiation, viewed as a localization of the damage variable within the model. For the purpose of current investigation various values of w 1 are tested that either initiate damage or not and in the former case the possibility of localization is studied.

. . Problem setup

e problem setup is similar to what was described in Sec.( . . ) with some modi cations adapted to the current context only which are stated below.

e reference initial con guration Ω as shown in the Fig. ( . ) with the boundary ∂Ω now represents a porous solid instead of the porous skeleton as done earlier. At time t = 0 it is assumed that the porous solid is completely intact, stress free and fully saturated giving, α 0 (x) = 0, u 0 (x) = 0, ε 0 (x) = 0, S 0 (x) = 1, p 0 (x) = 0.

( . )

e mechanical boundary conditions are exactly the same as stated in Eq.( . ) and damage is free to evolve at all the boundaries whenever the damage criterion is met according to Eq.( . ).

For the hydraulic problem, the lateral and bo om faces are set to be impermeable.

e loading is done by a constant imposed outward ux, -q f e 2 , on the top face thus inducing a drying e ect. While the drying ux measured using the discharge rates in S ( ) is found to be a function of time, we choose to make a simplifying assumption of constant ux. ese boundary conditions ∀t ≥ 0 read,

∇p t • (-e 1 ) (x 1 =0) = 0, ∇p t • e 1 (x 1 =L) = 0, - κ η w kw (p t )∇p t • (-e 2 ) (x 2 =0) = q f , ∇p t • e 2 (x 2 =H ) = 0. ( . )
. .

One-dimensional base solutions

It can be seen from the problem setup that the geometry, material properties and initial conditions are invariant along x 1 -direction. e loading and boundary conditions on the top and bo om faces are as well invariant along x 1 -direction. e boundary conditions on the lateral faces demand vanishing gradients of the solution along the x 1 -direction. So, even though the problem does not render itself for obtaining an easy exact solution due to its nonlinearity, one can a priori anticipate the existence of a particular class of solutions that are homogeneous in the x 1 -direction and are dependent only on x 2 . Such solutions are termed here as base solutions. In fact to obtain the base solutions one just needs to solve the problem in one-dimension along the x 2 -direction with appropriate boundary conditions, instead of the full two-dimensional problem posed earlier. is is the purpose of this section. Rigorous mathematical proof of existence and uniqueness of such base solutions for all times is out of the scope of the current work. Instead we exhibit numerically these solutions for t > 0.

e domain of the one-dimensional problem is de ned Ω = (0, H ) along the x 2 -coordinate, with boundary ∂ Ω = {(x 2 = 0) ∪ (x 2 = H )}. e initial and boundary conditions from Sec.( . . ) are adapted to this domain. While studying the in uence of depth H could be interesting, we choose a depth such that H . In view of the damage criterion, Eq.( . ), with the de nition of average pore-pressure, p * t (S t , α t ) = p t (S t , α t ) S t , the solution states can be classi ed into two: undamaged (α t = 0 ∀x 2 ∈ Ω) and damaged (∃x 2 ∈ Ω | α t > 0). Since the initial state at t = 0 is assumed to be that of a uniformly intact solid, damage would initiate if the following local damage criterion is an equality. ) where Eq.( . ) gives -p t = p c (S t , α t ) = a (α t ) p c (S t ) and Eq.( . ) gives U (α t ; S t ) = a (α t ) U (S t ). In Eq.( . ) the terms that drive the damage evolution are rst and second. Now for α t = 0, a (α t ) = -2 and for 0 < S t < 1, U (S t ) > 0 and p c (S t ) > 0. Since ϕ t ∈ P the rst driving term is negative whose magnitude increases with decreasing S t and conversely with increasing α t . With a similar reasoning also the second term is negative. e third term is a positive constant, w (α t ) = w 1 , acting as the threshold for damage initiation. Owing to the drying ux applied at x 2 = 0 the saturation degree is the lowest here. us one can anticipate that starting from a uniformly undamaged state the damage initiates at x 2 = 0 within a nite depth and propagates into the domain driven by desaturation. However, at initial times since a minimum amount of capillary energy needs to be stored before damage initiation according to Eq.( . ), the saturation degree everywhere in the domain could be such that the damage criterion is a strict inequality and damage does not initiate at all. So, one can envisage a nite time t d > 0 beyond which the damage criterion is an equality at the drying boundary and so damage initiates. All the solution states for t ≥ t d are called the damaged states and this t d depends on the intensity of the drying ux for a given material. Figs. ( . & . ) show the time evolution of the solution and functions of the solution, with H = 1.0m, q f = 10m.s -1 , w 1 = 1000N.m -3 . e material properties correspond to sand and the uid combination to air-water. It can be seen that initially damage does not evolve anywhere in the domain even though ux induced desaturation initiates at x 2 = 0 and progresses into the domain. During this stage of evolution the strain remains compressive and the porosity reduces from the initial value as expected within the unsaturated zone. However, once the damage criterion becomes an equality, damage initiates at x 2 = 0. With damage growing from 0, the saturation degree degrades close to 0 indicating invasion of air phase while the pore-water pressure only reduces slightly in magnitude.

ϕ t U (α t ; S t ) -p t S t ϕ t + w (α t ) ≥ 0 in Ω, ( .
e path followed by the solution, in the space (S t , p t ), at x 2 = 0 and within the domain close to this boundary are shown in Fig. ( . (a)). It is clear that the path followed during the initial damage evolution is that of a 'so ening' with respect to pore-water pressure. Also interesting to observe is that even though the drying ux at the boundary continues to drive the desaturation process, the damage stagnates at a maximum value at larger times and this is replicated also within the domain, thus creating a zone of uniform damage. is is due to the simultaneous reduction in magnitude of the rst driving term in Eq.( . ) with increasing damage so that a er a certain value of damage, the criterion does not become an equality anymore even though desaturation progresses. Fig. ( . ) shows the time evolution of damage for di erent thresholds w 1 . In all cases the damage initiates at x 2 = 0 as anticipated and then propagates into the domain supported by a nite depth, d t , that increases with time. Lower thresholds are characterized by earlier initiation of damage and a damage value closer to 1 at the drying boundary.

. .

Bifurcation from and instability of the fundamental branch e one-dimensional base solutions resolved in Sec.( . . ) are the x 1 -homogeneous solutions or the so-called fundamental states of the corresponding twodimensional problem. ese solutions are denoted by (ū t , φt , ᾱt ) at time t and, for consistency, the solution of the hydraulic problem by ( St , pt ). e branch of the evolution along which all the states remain x 1 -homogeneous is called the fundamental branch. Following the works of B M ( ); P . ( b); S . ( ) for standard gradient damage models, we investigate the possibility of bifurcation from the fundamental branch of evolution. A bifurcation in this sense is understood as the availability of admissible solution state/s other than the fundamental one. Whether or not the solution shi s from the fundamental branch depends on the loss of stability of the fundamental state itself and the characteristics of the perturbation that can cause such a shi . In the earlier mentioned works concerning standard gradient damage models, the possibility of bifurcations was associated to a loss of uniqueness criterion of the fundamental state. We follow a similar approach, by rst analyzing the loss of stability of the fundamental state and then the possibility of bifurcation.

Loss of stability

By construction the fundamental state, (ū t , φt , ᾱt ), satis es the rst-order stability, Eq.( .). e question is, if it satis es the full stability condition. To verify this one needs to analyze the condition Eq.( . ) that is obtained by the expanding the perturbed state of energy in the vicinity of the one associated to the fundamental state,

hE t (ū t , φt , ᾱt )( -ūt , φ -φt , β -ᾱt ) + h 2 2 E t (ū t , φt , ᾱt )( -ūt , φ -φt , β -ᾱt ) 2 + o(h 2 ) ≥ 0. ( . )
By virtue of the equilibrium, Eq.( . ), and zeroth-order balance, Eq.( . ), the FDD in Eq.( . ) reduces to, ) with the directions (ūt , φφt , βᾱt ) ∈ C × P × D d t , with P de ned in Eq.( . ) and ) where Ω d t = (0, L) × (0, d t ) is the damaged strip if damage initiates and Ω e t = Ω\Ω d t is the undamaged region.

E t (ū t , φt , ᾱt )( -ūt , φ -φt , β -ᾱt ) = ∫ Ω -π t φt + w ( ᾱt ) (β -ᾱt ) + w 1 2 ∇ ᾱt • ∇(β -ᾱt ) dx, ( .
D d t = { β ∈ H 1 (Ω) : 0 ≤ β < 1 in Ω d t , β = 0 in Ω e t }, ( .

Undamaged states:

For t < t d , ᾱt = 0 throughout Ω and the FDD is positive for all directions (ūt , φφt , β) if β 0. For directions such that β = ᾱt = 0, the FDD vanishes and the SDD, Eq.( . ), reduces to,

E t (ū t , φt , ᾱt )( -ūt , φ -φt , β -ᾱt ) 2 = ∫ Ω ε( -ūt ) . . C . . ε( -ūt ) + N bϵ( -ūt ) -(φ -φt ) 2 dx . ( . )
which is positive. Hence, full stability of the fundamental state is guaranteed for t < t d .

Damaged states:

For t ≥ t d , the damage criterion is an equality within the damaged strip Ω d t = (0, L) × (0, d t ) and an inequality everywhere else, Ω e t = Ω\Ω d t . If we consider directions such that β ᾱt then the FDD is positive within Ω e t owing to the strict inequality of the local damage criterion. Whereas a vanishing FDD corresponds to directions such that β = ᾱt = 0 in Ω e t . In such directions stability of the fundamental state is guaranteed if, . ) and only if the above the inequality is not strict. Note that in the above the expression of SDD is given by Eq.( . ) and not by the reduced Eq.( . ) owing to the possibility of β ᾱt in Ω d t . D d t ⊂ D t , is the set of admissible damage elds given by, Noting, for the moment, that SDD in the damaged states is evaluated at the fundamental state, analysis of the condition in Eq.( .) is deferred to a later moment in view of studying the possibility of bifurcations.

E t (ū t , φt , ᾱt )( -ūt , φ -φt , β -ᾱt ) 2 > 0, ∀( , φ, β) ∈ C × P × D d t , (
D d t = {β ∈ H 1 (Ω) : ᾱt ≤ β < 1 in Ω d t , β = 0 in Ω e t }. (
Synopsis: While the undamaged states are fully stable, the damaged states are conditionally stable and this condition amounts to verifying the positive de niteness of the quadratic form E t (ū t , φt , ᾱt )(•) 2 when applied on directions belonging to C × P × D d t . Speci cally,

∀( , φ, β) ∈ C × P × D d t E t (ū t , φt , ᾱt )( , φ, β) 2 > 0 =⇒ stability E t (ū t , φt , ᾱt )( , φ, β) 2 < 0 =⇒ instability. ( . )

Possibility of bifurcations

As mentioned earlier, the notion of bifurcation from the fundamental branch is intertwined with the availability of admissible solution state/s other than the fundamental one. is amounts to studying the evolution problem at time t + η, with η > 0 small enough, knowing that the solution has followed the fundamental branch up until t, i.e., knowing the solution at t is (ū t , φt , ᾱt ).

Before proceeding further we assume that such an evolution is smooth enough that forward time derivatives of the solution components exist at time t and are de ned as,

u = lim η→0 u t +η -ūt η , ϕ = lim η→0 ϕ t +η -φt η , α = lim η→0 α t +η -ᾱt η . ( . )
Further a smooth growth of the damage zone is assumed within the time interval (t, t + η). More detail of such a hypothesis can be found in S . ( ). Now, if the solution continues to stay on the fundamental branch then it implies that ( u, ϕ, α) are as well x 1 -homogeneous giving, ( u, ϕ, α) = ( ū, φ, ᾱ). However, if bifurcation is to happen then there should be some other solution rate ( u, ϕ, α). To check this possibility we derive from the evolution problem at t + η, a rate problem that any solution rate needs to satisfy.

Imposing (u t +η , ϕ t +η , α t +η ) to satisfy the three evolution principles in Sec.( . . ) we get:

(a) Irreversibility of damage: Damage must be non-decreasing, i.e., α ≥ 0.

Consequently ( u, ϕ, α) ∈ C × P × D, with,

P = {φ ∈ H 1 (Ω)}, D = {β ∈ H 1 (Ω) : β ≥ 0 in Ω}. ( . )
(b) Stability: Directionally stability implies rst-order stability which at (u t +η , ϕ t +η , α t +η ) reads,

E t +η (u t +η , ϕ t +η , α t +η )( , φ, β) ≥ 0 ∀ ( , φ, β) ∈ C × P × D. ( . )
As analyzed earlier for loss of stability of the fundamental state at time t, for directions ( , φ, β) such that E t (ū t , φt , ᾱt )( , φ, β) > 0, by continuity, Eq.( . ) and consequently full stability of the state (u t +η , ϕ t +η , α t +η ) hold true.

Whereas for directions such that E t (ū t , φt , ᾱt )( , φ, β) = 0 this is not evident. In fact such directions correspond to β = 0 in Ω e t , i.e., β ∈ D d t . Dividing Eq.( . ) by η and passing to the limit when η tends to 0 gives the following condition on any ( u, ϕ, α)

∈ C × P × D d t , E t (ū t , φt , ᾱt ) ( u, ϕ, α), ( , φ, β) + E t (ū t , φt , ᾱt )( , φ, β) ≥ 0 ∀ ( , φ, β) ∈ C × P × D d t , ( . )
with , ( . )

D d t = {β ∈ H 1 (Ω) : β ≥ 0 in Ω d t , β = 0 in Ω e t }. (
Note the restriction of α to D d t as opposed to D done a few steps earlier when imposing irreversibility of damage. is is a consequence of the restriction imposed on the directions to be considered for E t (ū t , φt , ᾱt )( , φ, β ) = 0 to hold true.

A rst expansion of E t +η (u t +η , ϕ t +η , α t +η ) on the le -hand side of the above equation up to second-order leads to, ) where • denotes the norm of C × P × D. A second expansion in time up to second-order of the operators at t + η leads to,

E t +η (ū t , φt , ᾱt ) + E t +η (ū t , φt , ᾱt )(u t +η -ūt , ϕ t +η -φt , α t +η -ᾱt ) + 1 2 E t +η (ū t , φt , ᾱt )(u t +η -ūt , ϕ t +η -φt , α t +η -ᾱt ) 2 + o( (u t +η -ūt , ϕ t +η -φt , α t +η -ᾱt ) 2 ) = E t (ū t , φt , ᾱt ) - t +h ∫ t ∫ Ω ∂π s ∂S s ϕ s S s dx ds, ( .
E t (ū t , φt , ᾱt )(u t +η -ūt , ϕ t +η -φt , α t +η -ᾱt ) + η 2 E t (ū t , φt , ᾱt )( u, ϕ, α) + η 2 2 E t (ū t , φt , ᾱt )( u, ϕ, α) 2 + η E t (ū t , φt , ᾱt ) + η 2 2 E t (ū t , φt , ᾱt ) + o(η 2 ) = -η ∫ Ω ∂π t ∂S t φt S t dx - η 2 2 d dt ∫ Ω ∂π t ∂S t φt S t dx . ( . )
e following de nitions are obtained by di erentiating w.r.t time the energy balance wri en at time t,

E t (ū t , φt , ᾱt ) = - ∫ Ω ∂π t ∂S t φt S t dx, E t (ū t , φt , ᾱt ) = - d dt ∫ Ω ∂π t ∂S t φt S t dx .
( . )

e term E t (ū t , φt , ᾱt )(u t +ηūt , ϕ t +η -φt , α t +ηᾱt ) is o(η 2 ) due to the assumption of smooth growth of the damage zone and α ∈ D d t , see S

. ( ) for a detailed deduction. So, one obtains by diving Eq.( . ) by η 2 and passing to the limit when η tends to 0,

E t (ū t , φt , ᾱt )( u, ϕ, α) + E t (ū t , φt , ᾱt )( u, ϕ, α) 2 = 0. ( . )
e rate problem:

Subtracting Eq.( . ) from Eq.( . ) gives the following inequality that the rate at any time t > 0, ( u, ϕ, α) ∈ C × P × D d t , needs to satisfy,

E t (ū t , φt , ᾱt ) ( u, ϕ, α), ( -u, φ -ϕ, β -α) + E t (ū t , φt , ᾱt )( -u, φ -ϕ, β -α) ≥ 0 ∀ ( , φ, β) ∈ C × P × D d t . ( . )
To verify if any rate other than the x 1 -homogeneous rate, ( ū, φ, ᾱ), satis es the rate problem we proceed as follows. First choosing ( , φ, β) = ( ū, φ, ᾱ) in Eq.( . ) we get,

E t (ū t , φt , ᾱt ) ( u, ϕ, α), ( ū -u, φ -ϕ, ᾱ -α) + E t (ū t , φt , ᾱt )( ū -u, φ -ϕ, ᾱ -α) ≥ 0.
( . )

Whereas by choosing ( , φ, β) = ( u, ϕ, α) in Eq.( . ) that is already wri en with ( ū, φ, ᾱ) as the solution we get,

E t (ū t , φt , ᾱt ) ( ū, φ, ᾱ), ( u -ū, ϕ -φ, α -ᾱ) + E t (ū t , φt , ᾱt )( u -ū, ϕ -φ, α -ᾱ) ≥ 0. ( . )
Adding Eq.( . ) and Eq.( . ) we get the inequality, ) is condition holds true only if ( u, ϕ, α) = ( ū, φ, ᾱ), when E t (ū t , φt , ᾱt ) is positive de nite thus indicating a uniqueness criterion for the solution of the rate problem. Note that both ( u, ϕ, α) and ( ū, φ, ᾱ) come from the same space. Whereas, ( uū, ϕφ, αᾱ) ∈ C × P × D t where D t is the linear space associated to D d t at time t de ned as, ) Synopsis: e uniqueness of the x 1 -homogeneous rate solution of the rate problem and consequently the impossibility of bifurcation from the fundamental branch are ensured by the positive de niteness of the quadratic form E t (ū t , φt , ᾱt )(•) 2 when applied on directions belonging to C × P × D t . Speci cally,

E t (ū t , φt , ᾱt )( u -ū, ϕ -φ, α -ᾱ) 2 ≤ 0, ( .
D t = {β ∈ H 1 (Ω) : β = 0 in Ω e t }. ( .
∀( , φ, β) ∈ C × P × D t E t (ū t , φt , ᾱt )( , φ, β) 2 > 0 =⇒ no bifurcation E t (ū t , φt , ᾱt )( , φ, β) 2 ≤ 0 =⇒ bifurcation possible. ( . )
Remark: One can notice from the above results that both the loss of stability of the fundamental state and the possibility of bifurcation from it, amount to studying the positive de niteness of the same quadratic form, E t (ū t , φt , ᾱt )(•) 2 , when applied on directions belonging to di erent spaces: C × P × D d t and C × P × D t respectively. If t b and t s are two positive times at which respectively, bifurcation from and instability of the fundamental state rst occur, then one can prove as a consequence of P ⊂ P and D d t ⊂ D t that t b ≤ t s . is proof is considered to be out of scope of the current work. However, one can nd such a proof in S . ( ) employing the minimization of 'Rayleigh' ratios for both cases in a thermal shock problem which has a similar structure.

. . Characterization of bifurcations

Now we are in a position to study the possibility of bifurcations by assuming general forms of the directions on which the quadratic form is to be applied. To this end, let ( , φ, β) ∈ C × P × D t be a general direction. Its components can be decomposed into their respective Fourier modes with characteristic wave numbers k ∈ N as follows,

(x) = k k (x), k (x) = k 1 (x 2 ) sin kπx 1 L e 1 + k 2 (x 2 ) cos kπx 1 L e 2 , φ(x) = k φ k (x), φ k (x) = φ k 1 (x 2 ) cos kπx 1 L , β(x) = k β k (x), β k (x) = β k 1 (x 2 ) cos kπx 1 L .
( . )

While choosing the above decomposition, the boundary conditions that apply to the admissible perturbations to the fundamental state are a priori adapted to ( , φ, β). Speci cally, since the KKT conditions, Eq.( . ), demand that on any boundary where damage grows the spatial derivative normal to it vanishes. is applies to the lateral boundaries of the damaged strip and to the drying boundary, respectively {(x 1 = 0) ∪ (x 2 ∈ (0, d t ))} and {(x 1 = L) ∪ (x 2 ∈ (0, d t ))}. Whereas in Ω e t since damage is uniformly zero, its spatial derivative as well vanishes at its boundaries, {(x 1 = 0) ∪ (x 2 ∈ (d t , H ))} and {(x 1 = L) ∪ (x 2 ∈ (d t , H ))}. Consequently for the damage rate one gets, ∂ β/∂x 1 = 0 on x 1 = 0 and x 1 = L. Similarly using the boundary conditions on displacement imply that on x 1 = 0 and x 1 = L, 1 = 0 uniformly, implying ∂ 1 /∂x 2 = 0 and since the shear stress at this boundary vanishes we get ∂ 2 /∂x 1 = 0. In addition to the above, the x 1 independence of the fundamental state justi es the forms assumed for an β in Eq.( . ). Consequently the following de nitions for their gradients hold for each k, ( . )

ε( k ) = k 1 kπ L cos kπx 1 L e 1 ⊗
Owing to the classical constitutive relation for porosity, Eq.( . ), obtained by rearranging the zeroth-order balance law for variations in porosity, the above de nition of ε( k ) and the boundary conditions of the hydraulic problem for p t , Eq.( . ), one can justify the form assumed for φ in Eq.( . ). Exploiting the orthogonality among di erent Fourier modes allows to uncouple them and evaluate the functional E t (u t , ϕ t , α t )( , φ, β) 2 in Eq.( . ) at (ū t , φt , ᾱt ), in directions for each k separately. Following this and integrating once along x 1 gives,

E t (ū t , φt , ᾱt )( k , φ k , β k ) 2 = L 2 x 2 =H ∫ x 2 =0 λ k 1 kπ L + d k 2 dx 2 2 + 2µ       k 1 kπ L 2 + d k 2 dx 2 2      + µ d k 1 dx 2 -k 2 kπ L 2 + N b k 1 kπ L + d k 2 dx 2 -φ k 1 2 -φt π t (β k 1 ) 2 -2 π t φ k 1 β k 1 + w 1 2       β k 1 kπ L 2 + dβ k 1 dx 2 2      dx, ( . )
where functions denoted (•) are to be understood as evaluated at the fundamental state. Dependency on the wave number in the above expression can be understood as a parametrization and study of the positivity of quadratic form can be recast into a comparison with 0 the smallest eigenvalue of the eigenproblem, ) where µ k denotes the eigenvalues and ( k 1 , k 2 , φ k 1 , β k 1 ) the eigenvectors for each k.

E t (ū t , φt , ᾱt ; k) -µ k I ( k 1 , k 2 , φ k 1 , β k 1 ) 2 = 0, ∀k ∈ N, ( .
Accordingly the criterion for bifurcation, Eq.( . ), at each time t > 0 can be translated to,

∀k ∈ N ∀( k 1 , k 2 , φ k 1 , β k 1 ) ∈ H 1 ( Ω) × H 1 0 ( Ω) × H 1 ( Ω) × H 1 t ( Ω)
inf { (µ k )} > 0 =⇒ no bifurcation inf { (µ k )} ≤ 0 =⇒ bifurcation possible.

( . )

With Ω as de ned in Sec.( . . ), H 1 0 ( Ω) are the space of H 1 functions in Ω that vanish at x 2 = H and H 1 t ( Ω) are those vanishing in the undamaged region x 2 ∈ (d t , H ).

Indicative study of bifurcations

While the true possibility of bifurcations could be understood by the criterion Eq.( . ), we can investigate the behavior of the coe cients of the quadratic form E t (ū t , φt , ᾱt )( , φ, β) 2 which gives a qualitative indication. is is the purpose of the current section. By virtue of Eq.( . ) we can classify the various terms within the functional as follows, E t (ū t , φt , ᾱt )( , φ, β) 2 = A t (ū t , φt , ᾱt )( , φ, β) 2 + B t (ū t , φt , ᾱt )( , φ, β) 2 , ( . ) ( . )

We can conclude that the rst term, A t (ū t , φt , ᾱt )( , φ, β) 2 , is positive denite because N , w 1 are positive constants, the sti ness tensor, C, is positive de nite, φt ∈ P and π t = 2(-p c ( St ) St -U ( St )) < 0. us this term acts to prevent any bifurcations from occurring. e term B t (ū t , φt , ᾱt )( , φ, β) 2 on the other hand contributes to the possibility of bifurcations. Even if the sign of the cross term in ( φ β) could only be determined by solving the eigenproblem, its coe cient however is positive since π t = -2(1ᾱt )(-p c ( St ) St -U ( St )) > 0. For the sake of being conservative this cross term is considered negative in the following analysis.

In view of the above analysis, one can study the strength of B t by studying the strength of its coe cient as a function of time in order to understand the behavior of bifurcations, if they occur. It can be observed that this coe cient is function of the fundamental states that were resolved in Sec.( . . ). Fig. ( . ) shows the evolution of π t , with respect to the saturation degree for di erent values of the threshold w 1 . As the saturation degree reduces with time during the drying process it can be seen that this function deviates from its undamaged path as soon as damage is initiated and its magnitude tends to reduces as damage propagates into the domain. is indicates that the term B t reduces with time and one can infer that the tendency to bifurcate from the fundamental branch reduces as well. For higher values of w 1 however the magnitude of π t starts to pickup a er the initial reduction. However, the rate of this increase was found to be extremely slow due to the saturation degree being close to vanishing and the path followed is on a very steep degraded retention curve as seen in . One can make similar inferences for the loss of stability as well since, as remarked earlier, it also amounts to the study of a quadratic form involving the same operator, E t (ū t , φt , ᾱt ) according to Eq.( . ).

.

In this chapter a novel approach to modeling of capillary force driven fracturing phenomenon has been proposed, inspired by the experimental results mainly obtained by S S ( , a) and based in the gradient damage modeling following M . (

). e evolution problem for damage is posed starting from a total energy of the porous solid and following a variational approach assuming a minimization principle to hold true at each time. First-order stability conditions and necessary conditions are derived that form the governing equations that can be resolved numerically. Criterion for loss of stability of the base solutions and possibility of bifurcation from the fundamental branch are derived respectively from the second-order stability conditions and a rate problem. An indicative study reveals that, at the current state the model does not exhibit a clear bifurcation and similarly loss of stability of the damage base solution, which should induce damage localization, notwithstanding the behavior of π t which seems to indicate this possibility for relatively large values of w 1 .

While the solution of an eigenproblem could reveal more details, we claim that the lack of clear bifurcation is due to the simpli ed approach adopted in this paper. Speci cally, we did not account in a comprehensive manner for the degradation of the resistance to uid ow within the damaged zone and eventually along localized fracture planes. is could in principle be achieved via a suitable coupling between the intrinsic permeability and the damage parameter, for instance as proposed in M . ( ).

C e current thesis dealt with two types of instabilities: hydraulic and mechanical, that are intrinsic to partially saturated porous media such as geomaterials. It starts with a recount of the phenomenology of these two types of instabilities. en a concise description of two modeling approaches to partial saturation is done: a classical approach and a phase eld based approach. Following this, the thesis diverged into two parts each focusing on the one type of instability.

In Part II a one-dimensional analysis of the phase eld approach, assuming a rigid porous skeleton, revealed that moving fronts of invading uid phase can occur under suitable initial and boundary conditions. is is observed both in imbibition and drainage, thanks to the use of a nonconvex ux function. ese moving fronts are further resolved as traveling waves se ing the stage for a linear stability analysis.

e LSA and twodimensional simulations shed light on the hydraulic instabilities that can manifest themselves as uid ngering phenomena.

In Part III focused on mechanical instabilities in ne-grained geomaterials that result in formation of fractures typically with a complex pa ern. An extension of widely used gradient damage modeling of fracture to the context of partial saturation is done under the framework of a variational approach.

is modeling approach is tested in the case of a desiccating soil mass revealing its ability to describe fracture pa ern that is periodic as observed in representative experiments. While such a modeling perspective is the prevailing one, an investigation of a novel approach is done. is approach is inspired by experimental observations which have revealed that opening mode fractures of desaturating ne-grained soils is due to the invading uid lm rearranging the grains and not due to build-up of tensile stresses.

P is work leaves in its wake various possibilities, both modeling and experimental. e most natural one being the study of coupling between hydraulic and mechanical instabilities. e initiation of fractures induced by advancing drainage ngers is of particular interest in the context of CO 2 sequestration.

Concerning the phase eld partial saturation model, the rich structure of uid invasion involving coalescence and pinching needs cleverly designed experiments to validate. Also, the validation of results concerning the drainage front could be realized using quasi one-dimensional experiments that are scarce in existing literature. Coming to the capillarity induced damage model, while an indicative study of the bifurcation from fundamental branch is done, this modeling approach merits a full numerical resolution with material properties of various possible ne-grained soils.

Advanced numerical schemes are de nitely warranted both for the phase eld partial saturation model and the capillary gradient damage model, in view of accurate predictions.

A P P ENDICES *

In order to validate the robustness of the numerical scheme adopted to solve the one-dimensional problem of Secs.( . . & . . ) a mesh convergence analysis has been conducted for two test cases relative to the imbibition and the drainage problem. As expected reducing the element size implies a decrease of the -norm of the di erence between solutions, for the saturation degree and the generalized chemical potential, computed with successive mesh re nements, see Fig. ( . ). Data relative to the test cases are reported in the caption of the gure. (s e , 0, 0)

Starting from the system of equations Eq.( . ) rewri en as, S ξ (ξ ) = F (S(ξ )), ( . )

where S(ξ ) = [s, , w] T and F (S(ξ )) = [f ( ), (w), h(s, )] T is the associated short-hand notation for the right hand side of Eq.( . ). e system can then be expanded about an equilibrium S e (ξ ) = [s e , 0, 0] T as, S ξ (ξ ) = F (S e (ξ )) + DF (S e (ξ ))ϵ + D 2 F (S e (ξ )) . ) where DF (S e (ξ )) and D 2 F (S e (ξ )) represent the Jacobian and Hessian respectively, of F (S(ξ )) evaluated at the equilibrium. Subsequently ϵ = [s ϵ , ϵ , w ϵ ] T being a small perturbation in the vicinity of S e (ξ ) allows us to ignore O(ϵ 2 ) and higher order terms, leaving us the linear system of equations governing ϵ as, mis , numériquement, de modéliser des fronts de drainage, en plus de l'imbibition, sans employer de complexités supplémentaires dans le modèle. Une soigneuse analyse de stabilité linéaire (ASL) des solutions homogènes contre des perturbations arbitraires a été effectuée. La croissance dans le temps de l'amplitude des perturbations transversales imposées est comprise comme le phénomène de digitation. Les prédictions de l'ASL ont étés suivies par des simulations 2D mettant en lumière la capacité du modèle proposé de décrire le phénomène de digitation et la ségrégation des fluides.

ϵ 2 2! + . . . , (
Dans la deuxième partie, le déclenchement d'une fissure au sein d'un milieu poreux asséchant a été étudié. Dans l'esprit de la modélisation à gradient d'endommagement, un modèle de milieu poreux endommagé a d'abord été testé pour sa capacité à reproduire la formation de fissures périodiques telle qu'observée dans des expériences de laboratoire. A partir de ce cadre de la modélisation une généralisation de cette approche a été proposée en interprétant la présence d'une fissure dans un milieu poreux en train de sécher comme une perte de ses propriétés capillaires, permettant ainsi le passage d'un fluide non mouillant sous une pression capillaire évanescente. Ce paradigme est particulièrement intéressant dans la modélisation des sols à grains fins sans cohésion et non-consolidés, où la résistance à traction est négligeable et donc la fissuration induite en raison du développement des contraintes de traction n'est pas le phénomène dominant. Partant des principes qui caractérisent l'approche variationnelle utilisés pour construire le modèle, il a été montré que une dessiccation suffisamment forte peut amorcer l'endommagement de manière homogène sur la face en train de sécher tout en progressant dans le domaine dans le temps. L'apparition éventuelle de bifurcations de cette solution de base a été analysée, toujours dans le cadre de l'ASL.

Ce travail ouvre la voie à plusieurs pistes de recherche, la plus naturelle étant l'étude du couplage entre les instabilités présentées ci-dessus. En dehors de cela, des simulations bidimensionnelles de l'écoulement par digitation ont montré que le modèle prédit des caractéristiques supplémentaires d'écoulement instable, telles que le pincement et la coalescence de la phase de mouillage, qui doivent être étudiées à l'aide d'un campagne expérimentale soigneusement conçue. L'initiation de l'endommagement induit par l'évolution des doigts de drainage présente également un intérêt particulier dans le cadre des applications mentionnées précédemment. Enfin, des techniques numériques avancées peuvent être recherchées pour la résolution des problèmes ci-dessus dans le but de fournir des solutions précises plus efficacement.

Title: Study of hydro-mechanical instabilities in geomaterials

Keywords: partial saturation; geomaterials; phase field modeling; non-uniform fluid; gradient damage; Abstract: Hydraulic and mechanical instabilities in geomaterials refer to a variety of non-linear phenomena that can be triggered by heterogeneities inherent to such materials. Hydraulic instabilities in partially saturated conditions can manifest themselves as heterogeneous fluid invasion causing 'fingering' phenomenon. Mechanical instabilities on the other hand can present themselves as strain localizations and/or fractures. These instabilities and their associated coupling pose a major obstacle for a myriad of applications involving geomaterials like Carbon dioxide (CO2) sequestration, rain water infiltration and also for reliable predictions such as for contaminant flow in ground waters. Existing classical models do not resolve this behavior due to their lack of pattern-forming ingredients in their formulation and thus being stable against perturbations. The essence of current thesis work is to propose and investigate modeling techniques that allow to describe these instabilities in a simple and natural manner. The constitutive approach adopted in this thesis is that of micro-structured continua, in particular that of enhanced continua with a constitutive law depending on the gradient of so-called phase field variables.

The first part of the work deals with description of a fluid-fluid front that has been modeled as a diffused interface. This has been done by interpreting the presence of two fluids within the pore space as a single non-uniform fluid and the degree of saturation of one of the fluids as the corresponding phase field. While classical one-to-one relation between capillary pressure and saturation degree is expected to describe the retention properties of the porous network due to its texture, an enhanced relation is obtained by prescribing a chemical potential in the spirit of Cahn-Hilliard type modeling of multi-phase fluids. This enhanced relation together with a non-local energy contribution provides the ingredients required to describe fluid-fingering and non-local phenomena such as pinching and coalescence of fluids due to variations in capillary forces. Using asymptotic matching techniques it has been shown analytically, in a onedimensional setting, that the particular choice of energy contributions employed in the above framework allows to replicate 'overshoot' type non-monotonic saturation profiles during constant-rate infiltration into relatively dry medium. This has been found to match qualitatively with experimental observation. Further, a slight non-convexity introduced into the flux function has been shown nu-merically to allow the modeling of drainage fronts, besides imbibition, without employing any additional complexities into the model. A careful linear stability analysis (LSA) of the homogeneous solutions against arbitrary perturbations has been done. The growth in time of the amplitude of imposed transverse perturbations is understood as the fingering phenomenon. The predictions of the LSA have been followed-up by two-dimensional simulations portraying the ability of the proposed model to describe fluid fingering and segregation.

In the second part the triggering of a fracture within a drying porous medium has been studied. A prevailing modeling perspective, in the spirit of gradient damage modeling, has been first tested for its ability to replicate periodic fracture formation as observed in representative experiments. Further an investigation is done, starting from the gradient damage modeling framework, but interpreting the presence of a fracture within a drying porous material as a loss of its capillary properties, thus allowing the passage of non-wetting fluid under vanishing capillary pressure. This paradigm is of particular interest in modeling cohesion-less and unconsolidated fine-grained soils, where the resistance against tensile loading is negligible and thus fracturing induced due to development of tensile stresses is not the prevailing phenomenon. Starting from the principles of variational approach used to construct the model, it has been shown that for sufficiently strong desiccation, damage initiates homogeneously on the drying face while progressing into the body with time. The possible occurrence of bifurcations of this base solution has been analyzed, again in the framework of LSA.

This works sets the stage for various possibilities, the most natural one being the study of coupling between the above mentioned instabilities. Apart from that, two-dimensional simulations of fluid-fingering have shown that the model predicts additional features of unstable flow, such as pinching and coalescence of the wetting phase, which need to be investigated using carefully designed experiments. Initiation of damage induced due to evolving drainage finger is also of particular interest in the context of earlier mentioned applications. Lastly, advanced numerical techniques can be sought after for resolution of the above problems with an intent to provide accurate solutions more efficiently.

  phasic ow in porous media . . . . . . . . . . . . Fracturing in ne-grained soils . . . . . . . . . . . . . . . . . . Synopsis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Partially saturated poromechanics & a non-uniform uid . Classical approach to partial saturation . . . . . . . . . . . . . A phase eld approach to partial saturation . . . . . . . . .

Figure . :

 . Figure . : Schematic of geological storage of CO 2 . Inset: Schematic of core sample showing trapped CO 2 in yellow and the sealing cap-rock. Image provided by Global Carbon Capture and Storage Institute Ltd.

Figure

  Figure . : Dye tracer experiments (S ., ) exhibiting solute movement along preferential ow paths in relicit charcoal hearth soils.

Figure . :

 . Figure . : Phase diagrams from L () depicting qualitatively the zones of stable and unstable porous media ow. Irrespective of the we ing properties of the invading and defending uids these extremes of ow morphology were demonstrated experimentally to exist at similar extremes both for imbibition and drainage.

Figure .

  Figure . : e celebrated "Sa mann-Taylor"ngering observed at the interface between air (top) and glycerine (bo om), showing the development of instabilities and the inhibiting e ect that the fastest growing nger has on its neighbors (S T , ). In the rst image, the velocity of the interface was V = 0.1cm.s -1 and the critical wave length is c r = 1.2cm. e average wave length of the disturbance observed at initiation was 2.2cm which is in the expected range. Note that this experiment was done with the motion of the interface along the direction of gravity.

Figure . :

 . Figure . : Drainage experiments by M´ . ( ) where air (in white) displaced glycerol-water mixture (in gray) saturating a synthetic two-dimensional rigid porous medium. e drainage is simulated by constant rate with drawl from the bo om of the sample. e gravity component was tuned to study its stabilizing/destabilizing e ect, by tilting the sample from an inclined position ≈ 56 deg with respect to the ground (le -most), to a completely horizontal position parallel to the ground (right-most). disturbances of wavelengths larger than which destabilize the interface. b above is the gap between the plates of an analogous Hele-Shaw cell. e above analysis, complimented with experimental characterization (Fig.( . )) paved way to some pioneering works in the context of porous media ow. C . () proposed the existence of a macroscopic surface tension, akin to T in Eq.( . ), within porous media ow and thus having the same e ect as in Hele-Shaw ow. Later P F ( ) extended such result to three dimensions. Remark: While Hele-Shaw ow is mathematically akin to two-dimensional ow in porous media, there exist striking di erences especially in the transition region from one uid to another and in the morphology of the ngering instabilities (C W , ; S ., ). In porous media, the transition region from the part of the domain fully saturated with one uid, to the one which is fully saturated with the second uid can be described as a mixing region where both uids partially saturate the pores. On the other hand this transition is sharp in immiscible Hele-Shaw ow. is di erence is primarily due to existence of several internal length scales in the case of porous media (S ., ; H , ). erefore, the viscous ngering instabilities observed in porous media ow have a fractal structure rather than smooth formations as in the case of Hele-Shaw ow.

Figure . :

 . Figure . : Evolution of moisture content during in ltration into dry / silica sand measured by the intensity of light transmitted (G ., a). e light transmission through sands increases with an increase in the moisture content. Higher moisture content areas are thus visualized in the range of yellow to red.

Figure . :

 . Figure . : Results of saturation pro le versus in ltration ux, q [cm.m -1 ], in quasi one-dimensional experiments by (D C , ). e saturation pro les were obtained by light transmission through sands conned within tubes with diameter smaller than their typical nger sizes during in ltration.

Figure

  Figure . : Consolidated ger width versus in ltration ux results by D C ( ) of low (Y H , ) and high (G ., b) ux ranges with comparison of analytical nger width predictions by C . ( ); P H ( ).

Figure

  Figure . : Two-dimensional in ltration experiments by G . ( b) into initially dry / size fraction white silica sand.e experiments were conducted at increasing ratio, R s , of ux to saturated conductivity ux rates. R s = (0.012, 0.041, 0.088, 0.32, 0.82) to be read from le to right in the gure.

Figure

  Figure . : Schematic from S S ( b) of particle-level forces at play in the vicinity of a defect in the case of (a,b,c). miscible uid invasion and(d,e,f). immiscible uid invasion into an a priori saturated porous medium.

Figure . :

 . Figure . : Examples of desiccation induced fractures in di erent soils. Evenly spaced fractures in relatively shallow samples that are longer than they are wide in experiments with (a) clay (N K , ) (b) silt (P ., ). Polygonal fracture networks observed (P ., ) at the end of the drying experiment with square shaped silt samples of di erent depths: (c) 4mm, (d) 8mm.

Figure

  Figure . : (a) Experimental setup in S S ( ) for immiscible uid invasion. (b) Fractures originating at sub-millimeter "craters" on the surface of the sediment (immiscible uid invasion). (c) Fractures originating at the edges of an initial bore-shaped defect (miscible uid invasion). In both cases the fractures propagate vertically into the sample and laterally outwards.

Figure

  Figure . : Two-dimensional plane strain nite element simulations by S S ( a) of a representative Gaussian-shaped defect of depth 1mm within a rectangular domain (6 × 6)mm with an initial homogeneous void ratio of e(x, 0) = 3.92. e boundaries at the bottom and right-hand side are setup as no friction surfaces and the loading was a uniform suction on the whole domain increasing with time. (a) Void ratio, e(x, t ), when u c (t ) = 100KPa. (b) Plots depicting trends of the normal consolidation (NCL) and the air-entry (AEL) lines. e tip of the defect in the simulations is shown to deviate from NCL path thus reaching the AEL rst.

Figure

  Figure . : Pore-scale two-dimensional DEM simulations by J J ( ) with grain packs of 300 grains with size range [r min , 2r min ].e lateral boundaries were xed and the grain pack is compacted vertically under constant porepressure until vertical e ective stress of 3Mpa is reached. en an incrementally increasing gas pressure is applied at the bottom central pore. Maroon lines depict the compression at the inter-granular contacts, the blue dots are the pores fully occupied by gas and the pink lines are where the intergranular cohesive bonds are broken. (a) Classical capillary invasion for r min = 50µm at capillary pressure of 6KPa. (b) Fracture opening for r min = 0.1µm at capillary pressure of 2.55MPa.

Figure

  Figure . : (a,b,c) Experimental results from H . ( ) of invasion morphologies for varying injection ux (q), glass bead diameter (d ), and con ning weight (w ). e legend reads CF: capillary ngering, VF: viscous ngering, VF/CF: transition between the two former regimes and FR: fracturing. (d) Schematic of the experimental setup.

Figure . :

 . Figure . : Graphs of (a) interfacial energy and (b) capillary pressure functions given by Eqs.( . & . ) respectively, with S r es w = 0.1567, π 0 = 2840.91 Pa and the van Genuchten model parameter m = 0.685 typical to silica sand saturated by air-water mixture.

Figure . :

 . Figure . : Graphs of the relative permeability functions in Eqs.( . & . ) with the van Genuchten model parameter m = 0.685 typical to silica sand saturated by airwater mixture.

  Figure . :e double-well structure of the energy densities in (a) means that the pure uid phases can coexist, since they are the isopotential, (b), minima. All functions plo ed are dimensionless according to the scheme presented in Chp.( )

Figure

  Figure . : (a) e e ective energy density is no more a symmetric double-well function. (b) e pure uid phases are no more the isopotential minima of the e ective energy density. All functions plo ed are dimensionless according to the scheme presented in Chp.( )

  Figure.-Relative permeability functions K(S r ), Eq.( . ), for di erent values of van Genuchten parameters (a, b, c). K C :(1.175, 0.85, 2.0), K N C :(8.0, 2.0, 0.5). In ection point is shown as a circle on curve K N C (S r ).

  . ) where a, b and c are real constants obtained through ing with experimental results. As can be seen in Fig.( . ), this functional form has the possibility to model a purely convex behavior, K C (S r ), and as well a convex-concave behavior, K N C (S r ). In what follows, in order to account for consolidation of the porous skeleton (B , ) and to demonstrate in the context of the current phase eld model the possibility of modeling air-water displacements representing both imbibition and drainage, we employ the convex-concave function, K(S r ) = K N C (S r ), with a = 8.0, b = 2.0 and c = 0.5.

Figure . :

 . Figure . : Schematic of the imbibition simulations showing boundary conditions applied on the one-dimensional domain. e arrow indicates the direction of motion of the self-similar interface. e form of natural boundary condition V f (x = 0) is given in Eq.( . ).

  . ) is imposed. is has an e ect of perturbing the saturation degree at the le boundary towards S -> S 0 r , inducing imbibition, while the normal derivative of µ pf vanishes. e boundary at the right, x = 1, is drained with a Dirichlet boundary condition on µ such that, µ (x =1) = µ e (S 0 r ) + P(x = 1). ( . ) e results of simulations for various values of S -and λ = ρ w are shown in Fig.(

Figure . :

 . Figure . : One dimensional imbibition solutions of the coupled system Eq.( . ), for λ = ρ w , with initial condition S 0 r = 0.20 and boundary conditions corresponding to a constant rate of injection Eq.( . ) with S -= 0.96, 0.85, 0.80, 0.75, ∆t = E-, ∆x = E-, (a) Saturation degree, S r , (b) regularized e ective chemical potential, µ p f . Solutions are shown at di erent time-steps and restricted spatial range for clarity. Out of this range, solutions are continuous and uniform extensions up to their respective boundary conditions.

Figure . :

 . Figure . : One dimensional imbibition solution S r of the coupled system Eq.( . ), for λ = ρ w , S 0 r = 0.20 and S -= 0.85; (a) Solution at various time steps, (b) Solutions in (a) transformed into the TW-coordinate ξ = xct with c ≈ 2.63 given by Eq.( . ). Shock solution, Eq.( . ), for S -= 0.85 and S + = S 0 r = 0.20 is shown with triangle markers. Inset in (b) shows the solution, now with the nite element nodes depicted as black dots, and the shock solution in the immediate vicinity of ξ = 0.

  x = 0, by imposing gradually decreasing Dirichlet boundary conditions on S r and µ until the intended air saturation, S d r and µ d = µ e (S d r ) are achieved within a nite time, t d , and then those boundary conditions are kept constant for all t > t d . See schematic Fig.( . ).

Figure . :

 . Figure . : Schematic of the drainage simulations showing boundary conditions applied on the one-dimensional domain for t > t d . e arrow indicates the direction of motion of the self-similar interface.

  = μ e (S r 0 ) + P (x=1) μ e (S r d ) + P (x=0) (x=0)

Figure . :

 . Figure . : One dimensional drainage solutions of the coupled system Eq.( . ), for λ = ρ w , with initial condition S 0 r = 0.99 and boundary conditions corresponding to drainage reaching S d r = 0.20, 0.25, 0.30 in t d = 100∆t , with ∆t = E-, ∆x = E-, (a) Saturation degree, S r , with the pair of intermediate states (S a , S b ) ≈ (0.847, 0.365), (b) regularized e ective chemical potential, µ p f . Solutions are shown at the same time-step. Insets focus on the oscillatory junction between the bo om rarefaction wave and the uniform solution state (a) S b ≈ 0.365, (b) µ e (S b ).

FigFigure

  Fig.( . ) shows the solutions for various values of S d r and λ = ρ w . We note that the solution at later times is composed of an expanding part connecting S d r to a uniform state, S b , which is then connected by a sharper transition to another uniform state S a . S a then connects to S 0 r through a second expanding part. It is interesting to observe that for all values of S d r chosen, the solution se les down to the same values of S a ≈ 0.847 and S b ≈ 0.365, see Fig.( . ).e solution, µ pf , has a similar structure with uniform states corresponding to µ e (S a ) and µ e (S b ) with a di used transition in between.

  Fig.( . (b)). In Sec.( . . ) the corresponding TWsolutions of Eq.( . ) will be looked which are also non-classical in the sense of shocks violating the entropy condition such that Eq.(4.15). Whereas under a transformation of variable ξ = x/t, the part connecting S d r to S b tends to collapse into the corresponding rarefaction solution, Eq.( . ), of the hyperbolic equation, Eq.( . ), with S -= S d r and S + = S b , see Fig.( . (c)).

Figure . :

 . Figure . : One dimensional drainage solution S r of the coupled system Eq.( . ), for λ = ρ w , S 0 r = 0.99 and S d r = 0.20; (a) Solution at various time steps, (b) Solutions in (a) transformed into the TW-coordinate ξ = xct with c ≈ 3.41 given by Eq.( . ). Shock solution, Eq.( . ), for S -= S b ≈ 0.365 and S + = S a ≈ 0.847 is shown with triangle markers. Inset in (b) shows the solution, now with the nite element nodes depicted as black dots, and the shock solution in the immediate vicinity of ξ = 0. (c) Solutions in (a) transformed into the coordinate ξ = x /t . Rarefaction solution, Eq.( . ), for S -= S d r and S + = S b ≈ 0.365 is shown in circle markers.

Figure . :

 . Figure . : One dimensional drainage solution S r of the coupled system Eq.( . ), for λ = 1500ρ w , S 0 r = 0.99 and S d r = 0.20; (a) Solution at various time steps, (b) Solutions in (a) transformed into the TWcoordinate ξ = xct with c ≈ 7.05 given by Eq.( . ). Shock solution, Eq.( . ), for S -= S b ≈ 0.791 and S + = S 0 r = 0.99 is shown with triangle markers. Inset in (b) shows the solution, now with the nite element nodes depicted as black dots, and the shock solution in the immediate vicinity of ξ = 0. (c) Solutions in (a) transformed into the coordinate ξ = x /t . Rarefaction solution, Eq.( . ), for S -= S d r and S + = S b ≈ 0.791 is shown in circle markers.

Figure

  Figure . : (a)TW-solutions of Eq.( . ) for λ = ρ w , with boundary conditions s + = 0.20, s -= 0.96, 0.85, 0.80, 0.75, 0.70. Solutions are shown in a restricted range of ξ for clarity. (b) the corresponding connections shown on the curve K (S r ).

Figure

  Figure . : TW-solutions shown as trajectories in the phase space de ned by of Eq.( . ) for λ = ρ w , with boundary conditions s + = 0.20, s -= 0.96, 0.85, 0.80, 0.75, 0.70; phase plane views (a) (s, s ξ ), (b) (s, s ξ ξ ).

Figure

  Figure . : Implicit outer solutions, Eq.( . ), shown as solid lines. Corresponding shi ed solutions, ξ (s 0 ) -ξ (S c 1 ) and ξ (s 0 ) -ξ (S c 2 ), shown as dot-dashed lines. e jump part within the weak solution Eq.( . ) shown as do ed connection between points (0, S c 1 ) and (0, S c 2 ) which are shown as crosses. Boundary conditions shown as dot-dashed horizontal reference lines for (a) (s -, s + ) = (0.96, 0.20), (b) (s -, s + ) = (0.85, 0.20). Solutions are shown in a restricted range of ξ for clarity.

Figure

  Figure . : TW-solutions of Eq.( . ) that represent imbibition, (a) for varying s -, with s + = 0.20 and λ = ρ w , (b) for varying s + , with s -= 0.85 and λ = ρ w . (c) for varying λ, with s -= 0.85 and s + = 0.20. All the solutions are shown in a restricted range of ξ for clarity.

  Fig.( . (c)) shows the imbibition solutions for various values of λ and for xed boundary conditions, s -= 0.85, s + = 0.20.

Figure . :

 . Figure . : Trajectories of qualitatively di erent solutions with lower and higher starting equilibrium states, (S 0 b , 0, 0), for (a) λ = ρ w , (c) λ = 150ρ w . Zero slope locations on the trajectories are shown as circles and locations where the trajectory crosses the target S a are shown as lled dots. Solid lines represent the solution trajectory at the last iteration starting at equilibrium state (S k b , 0, 0). Corresponding TW-solutions centered within the domain ξ ∈ [-0.05, 0.05], for (b) lower values of λ = 0.5ρ w , ρ w , 1.5ρ w for which (s -, s + ) = (S k b , S a ) and S k b ≈ 0.30639, 0.36345, 0.413 respectively with λ and for (d) higher values of λ = 150ρ w , 1500ρ w for which (s -, s + ) = (S k b , S 0 r ) and S k b ≈ 0.79653, 0.77827 respectively with λ. Solutions are shown in a restricted range of ξ , s and aspect ratio of insets adjusted for clarity.

  Fig.( . (a)). For lower values of S 0 b

  Figure . : TW-solutions shown in Fig.( . b,d) shown as connections on the curve K (S r ) for (a) lower values of λ = 0.5ρ w , ρ w , 1.5ρ w for which (s -, s + ) = (S k b , S a ) and for (b) higher values of λ = 150ρ w , 1500ρ w for which (s -, s + ) = (S k b , S 0 r ) and s e 3 ≈ 0.81641, 0.83232 respectively with λ. e locations of (s e 3 , K (s e 3 )) on the curve are shown as circle markers.

  Fig.( . (b)) shows the TW-solutions thus obtained for di erent values of λ and the lower equilibrium state (S k b , 0, 0) associated to the last iteration.

  Fig.( . (c)). However, for lower values of S 0 b , a er oscillation, the trajectory reaches a local minimum without traversing the target and then diverges to full saturation. And for higher values of S 0 b , the trajectory reaches the target S a with a negative slope and then diverges down towards residual saturation. With these new stopping criteria for integration we obtain TW-solutions shown in Fig.( . (d)) for di erent values of λ. ese drainage TW-solutions are shown as connections on the permeability curve K(S r ) in Fig.(

Figure . :

 . Figure . : Dispersion curves of imbibition solutions for λ = ρ w , with boundary conditions s + = 0.20, a range of s -∈ [0.36, 0.96]. Respective s -values are labeled on each curve. Peaks that represent maximal positive growth rate are shown as lled dots on the curves, with their respective wave number, k f , labeled. Arrows trace the path followed by k f for decreasing s -.

Figure . :

 . Figure . : Plots indicating the correlation between the injection velocity and the nonmonotonicities in the base solutions and their associated instability behavior. (a)Wavelength associated to the fastest growing wave number, λ k (k f ), and the corresponding growth rate, σ s (k f ), as functions of the magnitude of dimensionless injection velocity, K (s -). (b) Maximum value a ained at the tip of the imbibition front in the base solutions, s max and the intensity of the nonmonotonicity represented by (s max -s -), as functions of K (s -).

FigureFigure . :

 . Figure . : (a) Eigen functions, s, for k = 140, λ = ρ w , with boundary conditions s + = 0.20, s -= 0.41, 0.50, 0.65, 0.85.e corresponding base solutions, s 0 , shown with triangle markers, are also plo ed here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.

  Fig.( . ) shows the dispersion relation for base solutions representing

Figure . :

 . Figure . : Plot indicating the correlation between initial saturation degree, s + , in the imbibition solutions and their associated instability behavior. Wavelength associated to the fastest growing wave number, λ k (k f ) is plo ed as a function of s + ∈ [0.16, 0.60].

  Fig.( . ) shows the dispersion relation for base solutions in Fig.( .

  Fig.( . ) shows the normalized eigen functions, s, for k = 140 along with the corresponding base solutions. It can be inferred that for larger values of λ (accordingly larger C a ) the eigen function has a peak value which is

Figure . :

 . Figure . : Dispersion curves of imbibition solutions for λ = 0.5ρ w , ρ w , 1.5ρ w , with boundary conditions s -= 0.85, s + = 0.20. Peaks that represent maximal positive growth rate are shown as lled dots on the curves, with their respective wave number, k f , labeled.

Figure . :

 . Figure . : (a) Eigen functions, s, for k = 140, λ = 0.5ρ w , ρ w , 1.5ρ w , with boundary conditions s -= 0.85, s + = 0.20.e corresponding base solutions, s 0 , shown with triangle markers, are also plo ed here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.

  Fig.( . ) shows the dispersion relation for base solutions in Fig.( . (b) & (d)) with varying λ. As mentioned earlier the drainage base solutions represent two qualitatively di erent type of TW-solutions.

  Figure.: Dispersion curves of drainage solutions for λ = 0.5ρ w , ρ w , 1.5ρ w , 150ρ w . Peaks that represent maximal positive growth rate are shown as lled dots on the curves, with their respective wave number, k f , labeled.

Figure . :

 . Figure . : (a) Eigen functions, s, of the drainage solutions for k = 150, λ = 0.5ρ w , ρ w , 1.5ρ w .e corresponding base solutions, s 0 , shown with triangle markers, are also plo ed here for reference, (b) close-up view showing the structure of eigen functions in the vicinity of their peaks.

  Figure . : (a) e relation Eq.( . ) plo ed in the space of (S 0 , k, (σ )). Only the unstable region, (σ ) > 0, is shown. (b) Contour plot of the same relation projected on the plane (S 0 , k ). e unstable region satisfying Eq.( . ) is shaded grey. (c) Dispersion curves plo ed for uniform base solutions S 0 = 0.96, 0.85, 0.80, 0.75, 0.70.

  Fig.( . (c)) shows the dispersion relations in the space of (k, (σ )) for uniform saturation values, S 0 , observed as le boundary conditions, s -, in the base solutions corresponding to imbibition in Fig.( . ). For S 0 = 0.70, 0.75 there is a range of k for which (σ ) > 0.

.

  Imbibition fronts . . . . . . . . . . . . . . . . . . . . . . . . . . Drainage fronts . . . . . . . . . . . . . . . . . . . . . . . . . . Synopsis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure . :

 . Figure . : Schematic of the imbibition simulations showing boundary conditions imposed. e arrow indicates the direction of motion of the interface. e form of natural boundary condition V f (x = 0, ) is given in Eq.( . ).

  . ) is imposed. is has an e ect of inducing imbibition with a uniform saturation degree of S -while the normal derivative of µ pf vanishes. e boundary at the right, (x = 1, ), is drained with a Dirichlet boundary condition on µ such that µ (x =1, ) = µ e (S 0 r ) + P(x = 1, ). ( . ) e top, (x, = 0.3), and the bo om, (x, = 0), boundaries are impermeable. See schematic Fig.( . ) for a summary of these boundary conditions.

Figure

  Figure . : (a) Evolution of initial perturbation to an imbibition front for S -= 0.75, S 0 r = 0.20 on a domain (x, ) ∈ ([0, 1], [0, 0.3]), with time step size ∆t = . E-, = 0.2m and characteristic element size ∆x = . E-= /4. Sequence shown is for increasing time step, t s t ep . e non-linear evolution involves formation of ngering like instability akin to observations in in ltration experiments. (b) Solution plo ed along the length of the central nger in solid lines, S r (x, = 0.15), and within the gap between the ngers in dashed line, S r (x, = 0.23), at t s t ep = 6000.

Figure . :

 . Figure . : Dimensionless quantities corresponding to the imbibition front shown in Fig.( . ) at t s t ep = 6000. (a) regularized e ective chemical potential, µ p f , (b) e ective chemical potential, µ e and (d) non-local potential, µ N L . (c) ese quantities plo ed along the length of the central nger in the vicinity of its tip, (x ∈ [0.795, 0.84], = 0.15).e non-monotonicities in S r and µ p f are not apparent due to this restricted spatial view.

  Fig.( . ) shows such a perturbed initial condition for the case of S -= 0.75, S 0 r = 0.20 and its evolution for t > 0 under the above boundary conditions and λ = ρ w . e perturbations grow in time forming ngers within the domain that are similar to observations in in ltration experiments into initially dry sand (G ., b; S .,). is model indicates that these ensuing ngers are more saturated than the region behind. Also, the advancing tips of the ngers are of higher saturation than the region within the nger itself. ese observations are inline with the experimental saturation pro les of ngers detected using light transmission by G . ( a). e saturation pro le within the space between the ngers, see Fig.( . (b)), is signi cantly di erent from that within the nger itself and is also closer to that of the one-dimensional base solution. ese intricate comparisons of the saturation pro les are yet to be realized using meticulously designed experiments. e solution of µ pf along with its constituent local and non-local components is shown in Fig.( . ) at t st ep = 6000, along with a comparative plot within the vicinity of the tip of the central nger. It is clear that while the e ective chemical potential, µ e , is non-monotonic owing to the non-convex local uid energy, the non-local chemical potential, µ N L , compensates by following the convexity of S r , in accordance with Eq.( . ). Along the sides of the ngers, the fringe region, the solution indicates regions of higher µ e .is is a consequence of higher saturation degree within this fringe region compared to that of the initially dry region, S 0 r , see Fig.( . (a)). is fringe region as well coincides with the passage of the distinctive swell of the nger tip at an earlier time step. e linear stability analysis in Sec.( . ) indicates that as the injection velocity, V f , is increased from the intermediate range, the spacing between

Figure . :

 . Figure . : Evolution of initial perturbation to imbibition fronts for injection velocities corresponding to Eq.( . ) for S -∈ [0.70, 0.96] (to be read row-wise), S 0 r = 0.20 on a domain (x, ) ∈ ([0, 1], [0, 0.3]), for increasing time step, t s t ep .

Figure . :

 . Figure . : Schematic of the drainage simulations showing boundary conditions imposed. e arrow indicates the direction of motion of the interface.

  Owing to the phase eld parameter, S r , of the current model being representative of the saturation degree of the we ing uid, an injection ux of the form Eq.( . ) at the boundary would not be appropriate to induce drainage. Instead the imposed linear pressure, P, is allowed to drive the drainage front across the domain by choosing appropriate boundary conditions. See schematic Fig.( . ). For all t > 0, at the le boundary surface, (x = 0, ), Dirichlet boundary conditions are setup with S r (x = 0, ) = S - and µ (x =0, )

Figure

  Figure . : (a) Evolution of initial perturbation to a drainage front for S -= 0.20, S 0 r = 0.99 on a domain (x, ) ∈ ([0, 1], [0, 0.15]), with ∆t = . E-, = 0.2m and ∆x = . E-= /4. Sequence shown is for increasing time step, t s t ep . (b) Solution plo ed along the length of the second nger from top in solid lines, S r (x, = 0.1), and within the gap between the ngers in dashed line, S r (x, = 0.0775), at t s t ep = 250, 2000.

  Fig.( . (b)) t st ep = 250, the la er being consistent with the one-dimensional base solution for the imposed boundary conditions in the current problem, see Fig.( . (b)).

Figure . :

 . Figure . : Dimensionless quantities corresponding to the drainage front shown in Fig.( . ) at t s t ep = 2000. (a) regularized e ective chemical potential, µ p f , (b) e ective chemical potential, µ e and (d) non-local potential, µ N L . (c) ese quantities plotted along the length of the second nger from top in the vicinity of its tip, (x ∈ [0.57, 0.615], = 0.1).

Figure . :

 . Figure . : Evolution of initial perturbation to drainage fronts for gradient of imposed pressure distribution λ = ρ w , 1.5ρ w , S 0 r = 0.20 on a domain (x, ) ∈ ([0, 1], [0, 0.15]), at t s t ep = 2000.

Figure

  Figure . : (a) Evolution of initial perturbation to a drainage front for S -= 0.20, S 0 r = 0.99 on a domain (x, ) ∈ ([0, 1], [0, 0.15]), with ∆t = . E-, = 0.2m and ∆x = . E-= /4. Sequence shown is for increasing time step, t s t ep . (b) Solution, S r (x, = 0.1) and S r (x, = 0.0775),plo ed along the of domain in solid lines and lines respectively. e spatial gap between these two solutions represents the initial perturbation at t s t ep = 0 and its evolution for all t s t ep > 0.

Figure . :

 . Figure . : Reference con guration of the desiccation problem with mechanical boundary conditions depicted.

Figure . :

 . Figure . : Evolution of the damage variable within the full domain ( rst column) and along the top boundary, x 2 = 0m (second column). e time corresponding to each row are t = {0.1, 0.8, 0.9, 1, 1.2, 10}s.

Figure . :

 . Figure . : Spatial evolution of the solution along the depth of the domain at the location x 1 = 0.069m and at time t = 0.8s. (a) strains, (b) stresses, (c) pressure and saturation degree of the we ing uid, (d) damage variable.sizes are ∆x = 0.001m and ∆t = 0.1s. At the beginning damage evolves in a homogeneous manner along x 1 driven by desaturation. At t = 0.9s the solution starts to bifurcates with a periodic structure with almost 10 sinusoidal pro les along the length. As time passes damage grows within these sinusoidal pro les while the base of the pro les not showing any change.is leads to localization within these pro les with a so-called period doubling phenomenon observed in similar applications of gradient damage modeling (B ., ; S ., ; M ., ). In essence selective initiation and propagation of fractures occur resulting a complex fracture network in higher dimensions. ese localized fractures propagate further into the domain as the desaturation progresses.e periodicity and spacing is qualitatively similar to what was observed in experimental desiccation tests as noted in Sec.( . . ).When it comes to the driving force behind the damage growth in this problem setup, one needs to look at the x 1 -homogeneous solution structure prior to any bifurcation. Fig.( . ) shows the solution pro les along the depth of the domain at time t = 0.8s. While the pressure and the saturation pro les at the top boundary, x 2 = 0m, con rm to the boundary condition imposed, the solution within the pro le is only a ected up to a certain depth controlled by the uid ow problem and this depth grows in time. e compression induced due to this desaturation is apparent looking at the pro le of ϵ -, whereas ϵ + remains vanishing. ere is as well a signi cant contribution due to the deviatoric strain, ε D , which in fact helps to drive the damage in accordance with Eq.( . ). is can be seen in Fig.( . (d)) where damage

Fig

Figure

  Figure . : Spatial evolution of the components of the volumetric uid ow vector V at time t = 10s, (a,b) within the domain and (c,d) along the length of the domain at the top x 2 = 0.05m, behind the tip of the fractures, x 2 = 0.045m and ahead of the tip, x 2 = 0.04m. e center of the fractures at t = 10s are approximately depicted as vertical red lines across the plots (c) & (d).

Figure . :

 . Figure . : Conceptual schematic of a twodimensional drying porous network along with the development of an opening mode fracture network originating on the drying face. Inset: zoomed-in view of the immediate vicinity of a conceptualized fracture

Figure . :

 . Figure . : Graphs of (a) capillary pressure and (b) interfacial energy functions given by Eqs.( . & . ) respectively, with S r es w = 0, π 0 = 2840.91 Pa and the van Genuchten model parameter m = 0.685 typical to silica sand saturated by air-water mixture.

Figure . :

 . Figure . : Evolution in Fig.( . ) shown within a restricted computational domain, Ω = (0, 0.05m), close to the drying boundary for clarity.

Figure . :

 . Figure . : Evolution of damage for values of w 1 (a) 2000N.m -3 (b) 1000N.m -3 , (c) 500N.m -3 (d) 100N.m -3 .

  . ) e directions (ūt , φφt , βᾱt ) are further denoted as ( , φ, β) ∈ C × P × D d t , with an abuse of notation w.r.t the notation introduced in Sec.( . . ).

  with,A t (ū t , φt , ᾱt )( , φ, β) 2 = ∫ Ω ε( ) . . C . . ε( ) + N (bϵ( ) -φ) 2 -φt π t β 2 + w 1 2 ∇ β • ∇ β dx, B t (ū t , φt , ᾱt )( , φ, β) 2 = -∫ Ω 2 π t φ β dx .

  Figure . : Evolution of the function π t with St at the boundary x 2 = 0m and two locations within the domain close to this boundary. e di erent values of w 1 used are (a) 2000N.m -3 , (b) 1000N.m -3 , (c) 500N.m -3 , (d) 100N.m -3 . In each plot the possible evolution if damage remains vanishing, π (S w , 0), is also plo ed for reference.

Fig

  Fig.( . (a)). One can make similar inferences for the loss of stability as well since, as remarked earlier, it also amounts to the study of a quadratic form involving the same operator, E t (ū t , φt , ᾱt ) according to Eq.( . ).

Figure . :

 . Figure . : Mesh convergence for one-dimensional simulations in Sec.( . ) with λ = ρ w , ∆x /∆t = 50. | |∆(•) | | 2 represents the -norm of the di erence between solutions computed with successive mesh re nements. (a) Imbibition simulations with S 0 r = 0.20, S -= 0.80, ∆x = . E-. E-. E-, E-, E-. (b) Drainage simulations with S 0 r = 0.99, S d r = 0.20, ∆x = . E-, E-, E-, E-, E-.

•

  s e , 0) = C a (ϕc -K (s e ))δ R ϕD N L K(s e ) ; h (s e , 0) = µ e (s e ) ϕD N L .( . )e characteristic polynomial, p(β), whose roots are the eigen values associated with the system Eq.(.) is thus given by,p(β) = det[DF (S e (ξ )) -βI] = β 3 -h (s e , 0)β -h s e ,0). tensor product between two vectors • • • inner product between two tensors involving a single contraction • . . • inner product between two tensors involving a double contraction • . . . • inner product between two tensors involving a triple contraction • × • Cartesian product of two function spaces O + or • + positive part of a tensor or an operator or a functional (as de ned when invoked) • -or • - negative part of a tensor or an operator or a functional (as de ned when invoked) det • determinant of a square matrix or a second-order tensor tr• trace of a square matrix or a second-order tensor max(•, •) maximum among two arguments min(•, •) minimum among two arguments | • | absolute value • norm as de ned when invoked D Cd is document is composed using L A T E X macros and compiled using pdfT E X engine. e typeface used is Linux Libertine. e layout is inspired and adapted from the tu e-book, classic-thesis and two-thirds classes.

Titre:

  Étude des instabilités hydro-mécaniques dans les géomatériaux Mot clés : saturation partielle ; géomatériaux ; modélisation de champ de phase ; fluide non-uniforme ; gradient d'endommagement ; Résumé : Les instabilités hydro-mécaniques dans les géomatériaux font référence à une variété de phénomènes non-linéaires qui peuvent être déclenchés par les hétérogénéités intrinsèques de ces matériaux. Les instabilités hydrauliques dans des conditions partiellement saturées peuvent se manifester comme invasion hétérogène d'un fluide 'par digitation'. Les instabilités mécaniques, d'autre part, peuvent se présenter sous la forme de localisations de déformations et/ou de fissures. Ces instabilités potentiellement couplées constituent un obstacle majeur pour une myriade d'applications impliquant les géomatériaux comme la séquestration du dioxyde de carbone (CO2), l'infiltration des eaux pluviales et aussi pour des prédictions fiables telles que le flux de contaminants dans les eaux souterraines. Les modèles classiques existants ne résolvent pas ce comportement en raison de leur manque d'éléments capable de modéliser la formation d'une microstructure dans leur formulation et sont donc stables contre toute perturbation. Cette thèse vise à proposer et étudier des techniques de modélisation permettant de décrire ces instabilités de manière simple et naturelle. L'approche en loi de comportement adoptée est celle des milieux continus à micro-structure, en particulier celle des milieux continus équipés d'une loi de comportement dépendant du gradient de variables de type champ de phase. La première partie du travail porte sur la description d'un front fluide-fluide qui a été modélisé comme une interface diffuse. Cela a été fait en interprétant la présence de deux fluides dans l'espace poreux comme un seul fluide nonuniforme et le degré de saturation de l'un des fluides comme un champ de phase. Alors qu'on s'attend à ce la relation bijective classique entre la pression capillaire et le degré de saturation les propriétés de rétention du réseau poreux en raison sa texture, une relation enrichie est obtenue en prescrivant un potentiel chimique dans l'esprit de la modélisation de type Cahn-Hilliard de fluides multiphasiques. Cette relation enrichie, associée à une contribution énergétique non-locale, fournit les ingrédients nécessaires pour décrire les écoulement par digitation et les phénomènes non-locaux tels que le pincement et la coalescence des fluides dus aux variations des forces capillaires. En utilisant des techniques de développement asymptotique, il a été démontré analytiquement, dans le cadre d'un modèle unidimensionnel, que le choix des contributions énergétiques utilisées permet de reproduire des profils de saturation non-monotones de type 'overshoot' lors d'une infiltration à vitesse constante dans un milieu relativement sec. Cela s'est avéré correspondre qualitativement à l'observation expérimentale. De plus, l'introduction d'une légère non-convexité dans la fonction de perméabilité relative a a per-

.

  Unstable bi-phasic ow in porous media . . . . . . . . . . . . . E ect of gravity . . . . . . . . . . . . . . . . . . . . . . Fracturing in ne-grained soils . . . . . . . . . . . . . . . . . . . Drainage induced fracturing . . . . . . . . . . . . . . . Synopsis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  ( . ), a technique novel to this eld. e leading tips of the generated ngers were clearly of higher moisture content. is observation and related numerical studies (N

					.,	; E
	G	,	) have lead to the hypothesis that this non-monotonicity is the
	cause rather than the consequence of the ensuing instability (E
	G	,	). is 'overshoot' of water content has been further evidenced
	experimentally to be present during two-dimensional (B		.,	;
	S		F	,	) and quasi one-dimensional (D C	,	;
	D C		.,	) experiments with transverse domain size smaller than
	the nger size expected.

  is analyzed in detail for it ability to predict uid ngering phenomenon in Part II.

		:
	.	Summary of the model -Governing equations . . . . . . . .
	.	Dimensional analysis & scaling . . . . . . . . . . . . . . . .
		H Y DRAULIC
		I N S TA BILI TIES

. Limit of vanishing di usion & the ux function . . . . . . . . PDE simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . Imbibition . . . . . . . . . . . . . . . . . . . . . . . . . . Drainage . . . . . . . . . . . . . . . . . . . . . . . . . . Traveling wave analysis . . . . . . . . . . . . . . . . . . . . . . . Imbibition . . . . . . . . . . . . . . . . . . . . . . . . . . Drainage . . . . . . . . . . . . . . . . . . . . . . . . . . Synopsis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table . :

 . Dimensionless numbers corresponding to the material properties and model parameters in Table.( . ) used through Part II, unless mentioned otherwise.

	Dimensionless number	Value
	C a

.

  Mathematical model . . . . . . . . . . . . . . . . . . . . . . .

	. .	Gradient damage modeling setup . . . . . . . . . . .
	. .	Transient hydraulic problem . . . . . . . . . . . . . .
	. .	Evolution problem . . . . . . . . . . . . . . . . . . .
	. .	First-order stability & necessary conditions . . . . .
	. .	Tension-Compression split . . . . . . . . . . . . . . .
	. .	Fracture toughness & related assumptions . . . . . .

. Numerical approximation & Algorithm . . . . . . . . . . . . . Two-dimensional desiccation problem . . . . . . . . . . . . . . . Problem setup . . . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . Synopsis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Algorithm . . . . . . . . . . . . . Two-dimensional desiccation problem . . . . . . . . . . . . . . . Problem setup . . . . . . . . . . . . . . . . . . . . . . . . One-dimensional base solutions . . . . . . . . . . . . . . Bifurcation from and instability of the fundamental branch . . . . . . . . . . . . . . . . . . . . . . . . . . . . Characterization of bifurcations . . . . . . . . . . . . . Synopsis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

			'	'
	.	Mathematical model . . . . . . . . . . . . . . . . . . . . . . .
		. .	A 'damaged' porous solid . . . . . . . . . . . . . . .
		. .	Transient hydraulic problem . . . . . . . . . . . . . .
		. .	Evolution problem . . . . . . . . . . . . . . . . . . .
		. .	First-order stability & necessary conditions . . . . .
	.	Numerical approximation &

Table .

 .

	: Material properties, model pa-
	rameters used through Chp.( ), unless men-
	tioned otherwise.
	property/ parameter	Fine Sand
	ϕ

  Figure . : Evolution of the one-dimensional base solution and functions of the solution within the full computational domain, Ω = (0, 1m).

	10	10 -5									0		10 -4
	8
											-2
	6
	4
											-4
	2
	0										-6
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
	1										0
	0.8										-1000
	0.6
											-2000
	0.4
	0.2										-3000
	0										-4000
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
	1										0.3
	0.8
											0.2998
	0.6
	0.4
											0.2996
	0.2
	0										0.2994
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
	1	10 -4									0	10 -4
	0.8
											-2
	0.6
	0.4
											-4
	0.2
	0										-6
	0		0.02	0.04		0.06	0.08		0.1	0.12		0		0.02	0.04		0.06	0.08		0.1		0.12
	1										0
	0.8										-1000
	0.6
											-2000
	0.4
	0.2										-3000
	0										-4000
	0		0.02	0.04		0.06	0.08		0.1	0.12		0		0.02	0.04		0.06	0.08		0.1		0.12
	1										0.3
	0.8
											0.2998
	0.6
	0.4
											0.2996
	0.2
	0										0.2994
	0		0.02	0.04		0.06	0.08		0.1	0.12		0		0.02	0.04		0.06	0.08		0.1		0.12

  Energy Balance: e energy balance at time t + η reads,E t +η (u t +η , ϕ t +η , α t +η ) = E t (ū t , φt , ᾱt ) -

			. )
	(c) t +h ∫ t ∫ Ω	∂π s ∂S s	ϕ s S s dx ds.

  ⊗ e 2 + e 2 ⊗ e 1 ),∇β k = -β k

		e 1 +	d k 2 dx 2	cos	kπx 1 L	e 2 ⊗ e 2
	+	1 2 (e 1 1 sin kπx 1 L d k 1 dx 2 -k 2 kπ L kπ L sin kπx 1 L e 1 + dβ k 1 dx 2 cos	kπx 1 L	e 2 .

In the current work we make a simplifying assumption that the G c remains constant and the above relations are used to rephrase the total energy as,

.

&

When it comes to numerical resolution, in standard gradient damage modeling the minimization problem for the non-convex total energy potential is recast into two separate minimizations problems, w.r.t u and α, both of which are convex when the other variable is considered a constant. A popular algorithm for this purpose is the alternating algorithm proposed by B . ( ). In the current context given that the total energy is formulated considering the hydraulic problem as providing instantaneous input through the resolution of π , one can see that the alternating algorithm can be adapted for the minimization of E(u, α) given by Eq.( . ), w.r.t u and α separately subject to the unilateral constraint of irreversibility of damage and appropriate boundary conditions.

Given the solution triplet (u n-1 , π n-1 , α n-1 ) of the solid and uid problems at time-step (n -1), Alg.( ) describes the alternating algorithm to obtain the solution at time-step n. At each alternating iteration i, the nonlinear coupled problem for solid displacement and equivalent pore pressure, (u i n , π i n ), is resolved using an iterative approach (herea er referred to as the porosplit iterations) as mentioned earlier in Sec.( . . ). Speci cally this coupled problem is recast into two sub-problems at each poro-split iteration k: a linear coupled problem for (u i .), in tandem with the transient uid ow equation, Eq.( . ), both with appropriate boundary conditions. e time derivative within Eq.( .) is discretized using the implicit Euler scheme of rst-order with time-step size ∆t.

is coupled problem can be viewed as a saddle-point problem (M ., ) which demands a certain stability to the nite element approximation, especially in the case of undrained situations with Biot coe cient close to unity, b ≈ 1, in which case the problem is closer to in-compressible elasticity. Even though we do not intend to encounter such a limit case, the stability issue can be resolved by adhering to the celebrated LBB condition (B F , ), which can be practically implemented by choosing a nite element interpolation with higher polynomial degree for u compared to p. While other choices are possible, we employ a Taylor-Hood element pair with quadratic and linear Lagrange nite elements for u and p respectively. In order to resolve the coupled problem a sparse LU decomposition routine available in the FEniCS suite (A ., ) is employed. Once p i n k is resolved (S i n k ) is obtained trivially through the algebraic relation Eq.( . ) at the corresponding nodes. At the end of each poro-split iteration, the equivalent pore-pressure, π i n k , at each node is obtained using Eq.( . ). ese iterations are repeated until convergence of the 2 -norm of the maximum di erence among the solutions of successive iterations, max{

Algorithm : Alternating algorithm for capillary damage model

Input:

e above modeling approach is now applied to study the desiccation of soils, inspired by various desiccation experiments that can be found in the literature, P . ( ); S S ( , a); S ( ). For numerical purposes the ux on the drying face/s is estimated (S , ) using the discharge rate that is calculated from experimental measurements of mass loss of water. In line with what was done in Sec.( . ), we consider here as well a plane-strain assumption.

e boundary value problem formed by the coupled system of equations Eqs.( . , . ) and the bound-constrained minimization w.r.t α, Eq.( . ), are resolved with appropriately de ned boundary condition and initial conditions as laid out further in Sec.( . . ) and following the algorithm described in Alg.( ). Material properties of the porous medium and parameters of the model chosen for the purpose of simulations are listed in Table .( . ), which are in the range typical of ne silica sands saturated with air-water mixture.