Micromechanical formulation of first‐ and second‐order works in unsaturated granular media

The current work provides a micromechanical formulation of the first‐ and second‐order works for unsaturated granular media. For this purpose, a Representative Elementary Volume (REV) of a triphasic medium is considered, which consists of a collection of particles with the pore space in between them occupied by wetting and non‐wetting fluids. By assigning a kinematics to the REV at the micro‐scale, the first‐ and second‐order works of internal forces associated with the solid phase, fluids, and the interfaces between them are derived. The formulations are presented in the most general case considering large deformations within both Eulerian and Lagrangian formalisms, while the special case of infinitesimal deformation and rotation is also investigated. With respect to classical formulations, the obtained expressions for the first‐ and second‐order works include additional contributions related to the deformation of the interface of the two fluids, and in case of the first‐order work, the change in geometry of the intersection line of the three phases. Next, we consider a special case of pendular regime at low wetting saturation, where the wetting fluid appears as a set of isolated liquid bridges formed between particle pairs. Such condition has been previously simulated via Discrete Element Method (DEM) numerical framework where the distributed capillary forces are replaced with a resultant capillary force that is added as a new contact force. We show here that using the resultant capillary force, the generated first‐ and second‐order works in the REV are generally different from the ones produced by the distributed capillary forces, unless under special conditions.

disastrous events if not predicted in advance. It is thus of vital importance to have an understanding of the instability mechanisms in unsaturated soils and to develop accurate failure prediction criteria for geotechnical infrastructure health monitoring. An abundance of research has been dedicated to this challenging problem and various failure criteria have been proposed for unsaturated soils, the validity of which is under debate.
There are currently three main widely accepted criteria for describing failure in soils 1 : (1) the plastic limit which corresponds to vanishing of the determinant of tangent elasto-plastic matrix, (2) the localization condition, 2 and (3) Hill's instability condition. 3 The notion of ultimate failure in soils was long believed to be described by the first criterion (shear localized failure), but it is now well-established that the failure can occur well before this limit in both dry and unsaturated conditions, 4,5 as also confirmed in experiments. 6 Indeed, the third criterion is now known as the general and necessary condition for detection of instabilities 7 and is triggered before the other two criteria in loading paths that start from an isotropic initial state. 1 This criterion corresponds to a vanishing of the second-order work as a result of vanishing of the determinant of the symmetric part of the tangent elasto-plastic matrix. When this condition is met, in at least one loading direction, the material becomes unstable with a bifurcation. In the condition where second-order work vanishes in the current loading direction, a diffuse failure pattern is observed with a sudden degradation in strength. 8,9 Within the small strain assumption and in the dry case, the macroscopic second-order work is written at the material point level as 2  = based on the stress increment and the strain increment . The instability criterion is satisfied if a stress increment exists such that 2  ≤ 0. Micromechanical formulations of the second-order work criterion have been obtained for dry granular materials by Nicot et al. 7,8 , among others, based on the sum of second-order works at contacts between pairs of particles. In the quasi-static condition, such criterion reads 2 where  is the set of all contact points, is the incremental contact force, and is the relative displacement at the contact point of two particles. However, as mentioned in the above paragraphs, it is important to extend the classical second-order work criterion to incorporate the coupled action of capillary and mechanical forces in possible loading scenarios leading to an instability. Among the most important developments on the topic is the study by Buscarnera and di Prisco, 10 where a macroscopic expression for the second-order work is obtained based on mixture theory and energy considerations. A fundamental simplifying assumption in their approach is the superposed continua where the three phases of the unsaturated porous medium are assumed to be smeared over the domain of Representative Elementary Volume (REV). It is however interesting to tackle this problem from a micromechanical viewpoint, so that micromechanical aspects such as fabric can also be incorporated into the second-order work criterion. The current work pursues such development and formulates the first-and second-order works based on information at the contact scale. Instead of the superposed continua assumption, we consider the actual microstructure of the unsaturated REV. The first-and second-order works are then calculated by assigning a kinematics to the REV. To the best of authors' knowledge, no micromechanical formulation of the second-order work criterion for unsaturated soils has been presented in the literature. It is through such micromechanical formulation that the true mechanisms contributing to the second-order work, as the necessary condition for instability, can be identified for any unsaturated soil.
The micromechanical formulation of the second-order work requires first deriving an expression for the first-order work. This is done in the current work in the most general case considering large deformations within both the Eulerian and Lagrangian frameworks. As a secondary objective of the current work, the macroscopic stress in an unsaturated granular material is defined based on micromechanical features with the aid of the developed first-order work expression. The defined macroscopic stress coincides with the one obtained by Duriez et al. 11 via homogenization (i.e., direct integration of stresses within the REV), which was verified against DEM simulations. Also, in the dry case, the obtained stress expression reduces to the one obtained by Chang and Kuhn. 12 In the current work, the unsaturated granular medium is treated as a discrete system of rigid particles with the space in between them inhabited by the wetting and non-wetting fluids. Both the wetting and non-wetting phases are assumed to be incompressible. Figure 1 shows a schematic picture of the REV considered for the micromechanical analysis in the current work. Also, Figure 2 shows a particle pair and different geometrical features of the three-phase system that are used in the mathematical derivations in this work. It should be mentioned that Figures 1 and 2 correspond to a special case where the wetting phase forms isolated liquid bridges between particle pairs in the pendular regime; however, the derivations in this work are independent of the saturation condition of the REV, except in Section 4 where aspects of DEM modeling of unsaturated condition in the pendular regime are discussed.
The obtained expressions for the first-and second-order works show the need for consideration of additional terms with respect to the widely used classical formulations and recently introduced expressions, see Buscarnera and di Prisco 10 and Houlsby 13 for instance. Specifically, the additional terms in the first-order work expression are related to the work of F I G U R E 1 The REV of unsaturated granular material shown in the special case of pendular regime. Transformation from the initial configuration to the deformed configuration is defined via the deformation gradient tensor . REV, representative elementary volume.
F I G U R E 2 Geometry of two particles and the liquid bridge formed at their contact point.
surface tension in changing the geometry of contact line of the three phases and on the stretching of the interface between the wetting and non-wetting fluids. In case of the second-order work, the additional term is related to the second-order work of suction increment on deforming the interface between the wetting and non-wetting fluids. These additional terms are missing in the current formulations of first-and second-order works in the literature.
The final contribution of the current work pertains to the aspects of numerical modeling of the unsaturated condition in the DEM framework. The DEM simulation of the unsaturated condition is currently possible for an idealized assembly of rigid particles in the pendular regime where the wetting phase forms isolated liquid bridges between particle pairs. 11 Duriez et al. 11 used equivalent capillary forces in DEM replacing the actual distributed internal capillary forces. The equivalent capillary forces were defined such that the volume average of stress in the REV remains the same as in the case with actual distributed internal capillary forces. It is shown in the current work that this equivalent capillary force is equal to the resultant capillary force on particles. In this work, and based on the obtained expressions of the first-and second-order works, we aim at investigating the energy implications of representing the capillary interactions with a single resultant capillary force in DEM simulations of the unsaturated granular media in the pendular regime. It will be shown that using the resultant capillary force, the first-and second-order works generated in the REV will be generally different from the case where the actual distributed nature of capillary forces are considered, except in special cases where the rotation of particles is neglected, and in the case of second-order work with an additional constraint of small deformations. In other words, the resultant capillary force currently implemented in DEM cannot reproduce the first-and second-order work aspects of the actual distributed capillary forces, unless under special conditions. Finally, it should be mentioned that the current work is restricted to the following conditions: (1) rigid particles, (2) incompressible wetting and non-wetting fluids, (3) static equilibrium, (4) exclusion of body forces and moments, and (5) exclusion of contact moments. Considering these restrictions, the defined macroscopic stress is based on the assumption that a Cauchy continuum is sufficient to describe the macroscopic behavior of the REV. Extension of the work to higher-order continuum or Cosserat continuum descriptions require consideration of contact moments, among others, which is beyond the scope of the current work.
Throughout this paper, we use boldface symbols (such as ) to represent vectors and tensors, while regular symbols (such as ) denote scalars. Regular symbols with subscript are used to represent the components of vectors and tensors. Subscripts are typeset in italics and repetition of indices follows Einstein's summation convention. For instance, represents the component of vector , and represents the component of tensor . Both boldface and regular symbols can have superscripts, which do not follow Einstein's convention. The left superscript "0" refers to the initial undeformed configuration. The right superscript usually refers to the material phase(s); in that case, they are typeset in text mode. Throughout this manuscript, "w," "n," "p," and "f" represent the wetting fluid, the non-wetting fluid, the pore phase, and the fluid (either wetting or non-wetting) phase, respectively. Moreover, denotes the particles phase. As an example, 0 is the normal to the particle phase in the initial configuration. In case the right superscript is not typeset in text mode, its meaning will be explained upon definition. Finally, we use the sign convention of Continuum Mechanics in our derivations, i.e. compressive stress is taken to be negative.

2
FIRST-ORDER WORK

Equilibrium equations in the current configuration
The ensemble of rigid particles in the REV is denoted by the set , while the set  represents all the contact points between particles. The subsets  and  of  denote the internal and boundary contact points, respectively. The static equilibrium equations of particle ∈  in the current configuration (Eulerian) is expressed as the force and moment equilibria of Cauchy tractions acting on particle's surface Γ , which has an outward normal : In Equation (1), the tractions from the stresses within the interfaces are also taken into account. Γ , is the interface of phases α and , , , is the contact line where phases α, , and meet, , refers to the force exerted on particle α at the contact point , , is the inward conormal to Γ , , , is the surface tension of the interface of phases α and (see Figure 2), is the position vector relative to particle's reference point, and ijk is the permutation tensor. The equilibrium of forces and moments for the interfaces of particle and the fluids leads to the following relation 11 : where f represents the pressure of fluids. Using Equation (2) in Equation (1) and after some manipulations, we get the final form of the Eulerian force and moment equilibria relations of particles in the unsaturated F I G U R E 3 The resultant capillary force on particles (Equation (4)) is the summation of the surface tension force on the contact line (left hand side), and the capillary pressure on the wetted surface of the particles (right hand side). condition: where = n − w is the capillary pressure (matric suction). It is worth mentioning that the last two terms in Equation (3a) are the resultant capillary force (as opposed to the actual distributed capillary force) and equal to the equivalent capillary force derived based on stress homogenization and applied to particles by Duriez et al. 11 to mimic the unsaturated condition in discrete element simulations of granular materials: Equation (4) shows that the resultant capillary force on particles consists of two parts. One pertains to the capillary pressure on the wetted surface of the particles, and the other is the surface tension force on the contact line ,w,n (see Figure 3).
The fluids are assumed to be at rest; thus, the divergence of the shear part of their stress is zero. As such, the balance of linear momentum for the fluids in the current configuration (Eulerian) can be written as It is worth mentioning that there is no need to consider the equilibrium of interfaces for writing the first-and secondorder works as they are implicitly considered when writing the equilibrium equation of particles.

Equilibrium equations in the initial configuration
We denote the current and initial position vectors by and (see Figure 1), respectively. The deformation gradient tensor is defined as = ∕ with determinant = det . We recall the Nanson's formula d = −1 0 d 0 , which relates the surface vector d in the current configuration to the corresponding surface vector 0 d 0 in the initial configuration (see Lai et al. 14 , Section 3.27). Moreover, we use the following relation: which has a similar proof to that for the Nanson's formula (see Appendix A). Finally, the Cauchy stress and the first Piola-Kirchhoff (PK) stress Π are related through the Piola relations Π = −1 and = 1 Π . Using the above relations, the Eulerian force and moment equilibrium equations of particles presented in Equation (3) take the following forms when written in the undeformed configuration (Lagrangian): where the following first PK stresses are defined: For the fluid phases, the first PK stresses are defined based on Piola's relations as from which we get Using Equation (11) in Equation (5), the equilibrium equation of the fluids are written as Upon expanding Equation (12)

General Lagrangian description
The virtual (or first-order) work principle states that, for a system in static equilibrium, the work done by the external forces on the system is equal to the work done by the internal forces, for all virtual movements of the system from the equilibrium condition. In the current work, for the unsaturated granular REV, rather than postulating the principle of virtual work, we derive it as an extension of the (strong) equilibrium equations of particles and fluids by writing their weak forms. Such weak form is obtained by multiplying the force and moment equilibrium relations by their relevant virtual displacements and rotations, and integrating over the volume of all phases.
In case of discrete particles, the integration over volume is replaced by summation over all particles. For the particles, as shown in Figure 4, their rigid movement is described by the displacement of their reference point (center of gravity) , plus a rotation about the reference point. As such, for any point of the particle with relative position vector (with respect to the particle's reference point, see Figure 4), the displacement vector can be written as F I G U R E 4 Displacement in rigid particles based on translation and rotation of their reference point .
For the fluids on the other hand, a continuum virtual displacement f is defined in the bulk of the fluid with the only requirement that it satisfies the boundary conditions, that is, it conforms to the movement of the boundary of fluids. Although such continuum displacement field is not defined uniquely in the bulk of the fluid, it is straightforward to prove that the average strain, and specifically the total volume change of the fluid is defined uniquely through such non-unique displacement field. Most importantly, at the interface with particles, the virtual displacement of the particle and that of the fluid are continuous in the normal direction to the interface, while they are discontinuous in the tangential direction. In other words, we assume that there is sliding at the interface of particles and fluids.
In this way, the first-order work (per unit initial volume of the REV 0 ) expression  of the unsaturated granular assembly can be written by multiplying the force and moment equilibrium equations of particles and fluids (Equations 7 and 13) with the respective virtual displacement and rotation fields, and summing the four obtained equations as where 0 Ω is the domain of phase α. Equation (15) states that the total virtual work developed in the unsaturated granular assembly due to virtual displacement and rotation fields is equal to zero. This is a weaker condition of equilibrium than Equations (7) and (13), since the local (strong) condition of equilibrium for each phase is released, and a global equilibrium condition is imposed by using the virtual displacement and rotation fields as weighting coefficients. Equation (15) can be re-stated by re-arranging the terms as Next, we set to separate the contributions of external and internal forces in the first-order expression (Equation 16). Following the developments of Chang and Kuhn, 12 since the external work of contact forces should be based on their action on boundary displacements , which is generally different than the boundary displacement of peripheral particles , we add a term with zero contribution to Equation (15) (first term below) to get Moreover, the last two terms on the right-hand side of Equation (17) can be expanded as where the increment of deformation gradient is defined as = ( )∕ , and 0 Γ rev, is the portion of the boundary of REV which is occupied by phase . Substituting Equation (18) in Equation (17) and after some manipulations, including multiplication by 1∕ , the Lagrangian virtual work expression in terms of the separated contributions of external and internal forces is written as where and In Equation (21), at the contact point of particles and , the contact force c is the force exerted from on , while at the boundary of REV it is equal to the boundary force exerted at contact point . The relative displacement of particles at their contact point c is accordingly defined as In Equation (21), the relative incremental displacement̃is defined with respect to the displacement increment of particle as̃=

Eulerian description of first-order work
The Eulerian description of the internal and external first-order works can be obtained by replacing the expressions of first PK stresses (Equations 8-10) in Equation (21) and recalling the relations In the above relations, the incremental displacement gradient is defined as = .
Equation (24) is the most general form of the Eulerian internal first-order work expression, which is used in Section 2.4 to define a macro-stress for the unsaturated granular material. It should also be mentioned that Equations (24) and (25) can alternatively be derived by directly writing the weak form of the equilibrium equations in the current configuration presented in Section 2.1. Nevertheless, the Lagrangian description of the first-order work is first presented in Section 2.3.1, since it is required for the derivation of the second-order work in Section 3. For further simplification of Equation (24), we use the following relations: 1 where p is the total pore volume and is the saturation ratio of the wetting fluid. Moreover, using the Young-Laplace equation ( = w,n w ), assuming that w,n is constant on the fluids' interface, and the following equality: we can write Equation (28) can be demonstrated using Stokes' theorem following the same procedure as done in Duriez et al. 11 for Equation (20) in their work. Additionally, the relative displacement of the fluids and particles in the normal direction to their interface is zero (̃= 0) (see the beginning of Section 2.3.1); hence, for the particle-non-wetting fluid interface, we have n = . This together with the incompressibility of particles let us write Finally, we define the average stretch increment (per unit area) of the fluid-fluid interface aŝ In Equation (32), p is the pore volume fraction or porosity. Also, the first term on the right-hand side denotes the virtual work produced at contacts. The contact stress is classically referred to as the effective stress 15 for saturated soils, a state variable that, in conjunction with other internal state variables, controls the mechanical state of the saturated granular assembly. However, the existence of a state variable that plays the equivalent role of effective stress in saturated conditions is currently under debate for unsaturated granular media; nevertheless, assuming its existence, the contact stress is one of the candidates for defining such an effective stress in modern formulations. 16,17 As such, the first term in Equation (32) can also be regarded as the work of the effective stress. The second term refers to the work by the mean (or equivalent) pressure of fluids in expanding the pore volume. The third term corresponds to the work done by suction during a change in saturation. The fourth term denotes the work of surface tension in stretching the interface between wetting and non-wetting fluids (see Figure 5A). Finally, the last term corresponds to the work of surface tension in changing the topology of the intersection line of the three phases (see Figure 5B). The last term depends on the wetting angle, the current saturation, and the change in saturation. With respect to classical formulations, most notably the one by Houlsby 13 (also repeated in Buscarnera and di Prisco 10 ), the above expression of the first-order work includes two additional terms (the last two terms). This will be explained in more details in Section 2.4.

Macro-stress definition
In this section, we use energy aspects related to the first-order work in unsaturated granular medium to obtain an expression for the macroscopic stress in such medium. Also, we compare this expression with the ones presented in the literature.
Obtaining an expression for the macroscopic stress of a heterogeneous REV generally requires the knowledge of localization relations that can link the heterogeneous displacement field inside the REV to macroscopic (representative) kinematics of the REV. Such relations are straightforward to develop for dry granular media based on the displacement and rotation of the reference point of particles as done in Chang and Kuhn. 12 However, in the unsaturated case where the kinematics of liquids and their interfaces are involved, this process is generally not possible, unless with oversimplifying assumptions.
In the first part of this section, we follow a simplistic approach to obtain an expression for the macroscopic stress of the unsaturated granular medium. The aim is only to validate the first-order work expression presented in Section 2.3.2 (Equation 32). We show in the next part that the same stress description can be obtained following a more physically sound approach based on the virtual work principle.
As mentioned in the above paragraphs, the displacement field inside the REV of the unsaturated granular medium is highly heterogeneous. It is not possible to describe such displacement field in an analytical form. As such, we use the mean field (affine) kinematics approximation in the form: where and Ω are the increments of the macroscopic strain and rotation, which are determined based on best fit of and in Equation (33) on the displacement and rotation of rigid particles' reference points and boundary contact points. Equation (33) smooths out the highly heterogeneous displacement and rotation fields inside the REV and replaces the discrete system of particles with an artificial continuous system. It is important to re-iterate that Equation (33) only holds for the particles' reference points. In other words, only the position vector of particles' reference points and boundary contact points can be used in this equation.

Simplistic approach
Given that the contact moments between the particles and also on the boundary are assumed to be zero, the current unsaturated granular medium can be described as a Cauchy continuum, and the internal first-order work of the REV can be written based on the macroscopic strain increment and a macroscopic stress , which is equal to the volume average of stress within the REV, 12 as We describe the displacement of different points of particles using the displacement and rotation of their reference points (Equation 14) which upon using Equation (33) becomes For the boundary contact points, we have = For the displacements related to the liquid phases, we use an over-simplistic assumption and extend Equation (33) in those phases. Using these conditions for displacements throughout the REV, Equations (24) and (25) take the following forms: . (38) In the above equations, the relative position vector̃is defined with respect to the position vector of the reference point of particle as̃= − . In Equation (37), the branch vector c is defined as where at an internal contact point, it is equal to the relative position vector of reference points of particles and , while at a boundary contact point, it is equal to the relative position vector of the boundary contact point with respect to the particle's reference point. Moreover, The internal and external first-order works in Equations (37) and (38) must be equal for arbitrary choice of and Ω , which leads to the following relations: and Using Equation (42) in Equation (37) and comparing the resulting expression with Equation (34), we can get the macroscopic (average) stress of the REV based on either boundary tractions or internal force and stresses of the unsaturated REV: It is possible to further simplify the right-hand side of Equation (43) by expanding its last term. We use the Young-Laplace equation ( = w,n w ), assume that w,n is constant on the fluids' interface, and use the following equality 11 : which results in The obtained macroscopic Cauchy (Eulerian) stress of the REV from the first-order work expression has the same form as the one obtained via volume averaging by Duriez et al. 11 Although such agreement was expected, and seems rather obvious due to: (1) the use of mean field approximation, (2) assuming a Cauchy continuum, and (3) the extra simplistic assumption for the displacement field throughout the REV, it indeed validates the assumptions made in deriving the firstorder work expression, including the important assumption of not treating the interfaces as separate phases and excluding their equilibrium in Sections 2.1 and 2.2.

A more physical approach
In this part, we provide a more physical derivation of the macroscopic stress for unsaturated granular medium with the aid of the virtual work principle. As mentioned in Section 2.4.1, the macroscopic stress considered in this work is equal to the volume average of the stress within the REV, which using the divergence theorem can be written in terms of boundary tractions as The boundary tractions for the current REV include boundary contact forces and liquid pressures. Hence, Equation (47) can be expanded as which has the same form as the first term on the right-hand side of Equation (43). Using the discontinuous divergence theorem, the last two terms in Equation (48) can be stated as which gives the macroscopic stress as Next, we write the virtual work expression considering only the particles. Removing the last two terms in Equation (17), which are related to fluid phases and transforming the obtained relation to the current configuration, we get from which we can get the equality between internal work, and external work due to contact forces (considering only the particles) as Using the localization relations described at the beginning of Section 2.4.1 in Equation (52) we get Comparing the right and left sides of Equation (53) (external and internal works considering particles only), the term related to Ω should be equal to zero. Thus, we get (54) Replacing Equation (54) with Equation (50), it is straightforward to arrive at the same expression for the macroscopic stress as in Equation (43).
The derivation in this section does not involve the over-simplistic assumptions used in Section 2.4.1 and leads to the same macroscopic stress as in Duriez et al. 11 Moreover, the procedure shown in this section is important, as it gives the corrected version of the one presented in Li 17 with corrected virtual work expression for particles, and the correct use of the divergence theorem for converting the integral over surface tractions into volume integrals.

Discussion on energy aspects
In this part, we calculate the stored energy at contacts per unit volume of REV, which is also known as the work input to the unsaturated granular medium. 10 Setting Equations (32) and (34) equal and knowing that because of particles' incompressibility: p ∕ = , we get The above expression for the stored energy at contacts can be directly compared with the ones presented by Houlsby 13 and repeated in Buscarnera and di Prisco 10 where the last two terms are missing. Specifically, Houlsby 13 assumed that the interface between the wetting and non-wetting phases moves with the solid phase (particles) and as a result, the contribution of the associated surface tension disappears from the expression of the work input. Gray et al., 18 on the other hand, presented a more enriched expression than the one by Houlsby 13 for the work input to unsaturated elastic porous media based on thermodynamically constrained averaging theory. In their formulation, the work associated with the surface tension of the fluids was considered. However, there is no trace of the developed work at the contact line of three phases in their expression.
A more important discussion is regarding the definition of free energy in the unsaturated granular system, which leads to identifying work-conjugate state variables. The free energy is the mechanical energy that is not dissipated and is readily available to be restored. The contact term in Equation (55) is indeed the total mechanical energy provided to the particle assembly. In general, part of this energy is dissipated in form of heat due to friction, while the remaining part has two components. One component is the mechanical energy that is locked into the granular microstructure, and the other part is the mechanical energy that can be restored and is regarded as free energy. Neglecting dissipative mechanisms and the locked (frozen) energy at particle contacts, we can take the stored energy at contacts (Equation 55) as the free energy of the unsaturated system. At fixed saturation, a change in configurational strain leads to change in the topology of the liquid bridges. As such, the last two terms in Equation (55) are not only functions of the change in saturation degree ( ), but also the macroscopic incremental strain ( ). Such dependencies are very complex and as a result, it is not possible to provide analytical expressions for the work-conjugates to the macroscopic strain and saturation degree. This micromechanical investigation sheds light on the fact that the definitions provided in the literature for these variables, see, Buscarnera and di Prisco 13 ; Houlsby 10 , are not supported by micromechanics, and can only be viewed as first-order approximations.

Lagrangian description of the second-order work
In this section, we derive the Lagrangian form of the second-order work expression of the unsaturated granular medium by simply differentiating the first-order work expression derived in Section 2.3.1 (Equations 19-21), which is recalled here Since the first-order expression in Equation (56) is written in the initial configuration, it is straightforward to differentiate it paying attention to the fact that the configuration is fixed: In deriving Equation (57), the following Jacobi's formula for the variation of is used: Also, 2 represents the second-order variation with time. Using the Gauss theorem and the fluids' equilibrium equation (Equation 13), we can expand the following two terms from Equation (57) Moreover, we have Using Equations (58)-(60) in Equation (57) and using the particles' force and moment equilibrium equations (Equations 7a and 7b) together with the fact that  int =  ext , we obtain the final form of the second-order work expression in Lagrangian description as

Eulerian description of the second-order work
We Moreover, the relative displacement of the fluids and particle in the normal direction to their interface is zero (̃= 0). Also the particle is assumed to be incompressible. Using the above relations and after some manipulations, we obtain the expression of the second-order work in the Eulerian formalism, which after separating the contributions of external and internal forces becomes The displacement corresponding to the change in the topology of the interface between fluids (last term in Equation 70).
If we adopt the infinitesimal deformation and rotation assumption, the external and internal second-order works simplify as Equation (69) for the second-order work of external forces has the same general form as the one considered in Buscarnera and di Prisco 10 (Equation 5), except that instead of the assumption of superposed continua, our relation includes the actual incremental external forces at different parts of the boundary of the REV. In Equation (70), the first term on the right-hand side denotes the second-order work related to the solid skeleton that is developed at particles' contacts. The second term refers to the second-order work by the mean (or equivalent) pressure increment of fluids in expanding the pore volume. The third term corresponds to the second-order work done by suction increment during a change in saturation. The fourth term denotes the second-order work of suction increment in changing the topology of the interface between fluids (see Figure 6). With respect to classical formulations, most notably the one by Buscarnera and di Prisco, 10 the developed expression of the second-order work in Equation (70) includes one additional term (the last term) corresponding to the topology of the interface between fluids. This will be explained in more detail in the following.
Next, we re-write the developed second-order work (per unit volume) pertaining to the solid skeleton. The internal second-order work of the REV can be written alternatively based on a representative macroscopic strain increment and a macroscopic stress increment as Setting Equations (70) and (71) equal and knowing that because of particles' incompressibility: p ∕ = , we get Equation (72) can be directly compared with the expression of the second-order work proposed by Buscarnera and di Prisco. 10 . Apart from the factor 1∕2 considered in their derivation, the last term in our relation is absent from their result.

ASPECTS OF DISCRETE ELEMENT MODELING
In this section, we present aspects of simulating the unsaturated condition within the Discrete Element Method (DEM) framework. To this end, we consider the low wetting saturation condition (pendular regime) where the wetting phase forms isolated liquid bridges at the contact point of particles (see Figure 2). The neighboring particles sustain contact force , which is determined based on the assumed contact laws, that is, linear elastic or Coulombian friction, and is a function of the relative displacement of the particle pairs at their contact point. Scholtés et al. 19 and Duriez et al. 11 proposed extensions to DEM formulation for simulating the pendular regime for spherical particles within DEM. In their implementation, all the internal distributed forces of the unsaturated medium are described through interaction forces between particles. Specifically, the distributed capillary forces are lumped into an equivalent capillary force , which is derived based on homogenized stress of the REV and is equal to the resultant capillary force on particles derived in Section 2.1, and added to the contact force . Therefore, the wetting and non-wetting phases are not explicitly modeled in the DEM and their mechanical effects are only indirectly modeled through updating the equivalent capillary forces that are determined by computing the geometry of the liquid bridge. 11 As a result, this treatment of the unsaturated condition cannot model the pore fluid flow and the related transient effects. Nevertheless, Duriez et al. 11 showed, via numerical simulations, that such treatment of the unsaturated medium indeed gives an accurate estimate of the average stress in the REV.
In the current work, we investigate whether or not the equivalent capillary forces commonly implemented in DEM can fully represent all the effects of actual distributed capillary forces in terms of energy aspects. In other words, when the actual fluids are not explicitly modeled in the DEM, do the first-and second-order works generated in the REV based on the resultant capillary force match the actual ones based on distributed forces as obtained in the current work?

First-order work in DEM
The Eulerian first-order work expression was presented in Equation (32) The left-hand side terms have been collected to represent the first-order work associated with the commonly used, socalled net stress. It should be noted here that the current DEM model of unsaturated condition, which deals with resultant interaction forces, cannot describe the first-(or second-) order work in the biphasic condition, 11 as this condition produces zero net force on solid particles. It can only describe the work that is generated by the suction (the difference between nonwetting and wetting pressures) in the triphasic condition. This is the reason we look for the contribution of the resultant capillary force in the expression of the first-order work of the net stress (the right-hand side of Equation 73).
The second term on the right-hand side of Equation (73) can be expanded with the aid of Young-Laplace equation ( = w,n w ) and the equality = w as In the above relation, the last term on the right-hand side also includes the surface integrals related to the intersections of liquid bridges with the boundary of the REV, as those parts of the boundary are considered the surface of the external particles in contact with the internal peripheral particles. Moreover, it is important to mention that Equation (74) is only true for the special case of having isolated liquid bridges in the REV; otherwise, a surface integral term related to wetted surface in the boundary of the REV (Γ rev,w ) would appear on the right-hand side of this relation.
Using Equation (74) in Equation (73) and simplifying the resulting expression using Equation (28) gives . (75) Next, we look at the first-order work generated in the REV as a result of implementing the resultant capillary force in DEM: The resultant capillary force acting on the particle is presented in Equation (4). Looking at liquid bridge that connects two particles, since the liquid bridge is in equilibrium, such a capillary force is equal in magnitude on both sides of the liquid bridge. As a result, we can define a resultant capillary force for each contact, which we denoted in Equation (76) as , . We can turn the summation over contacts for the capillary forces into summation over particles, which after replacing the resultant capillary force from Equation (4) becomes Comparing Equations (75) and (77), we see that the calculated first-order work in DEM based on resultant capillary force does not match the first-order work obtained in the current work. This leads to an important conclusion that the simulation of unsaturated condition in DEM through equivalent capillary contact force has serious energy implications in the general case. Essentially, some part of the stored energy in the unsaturated granular medium will be lost following this approach. However, in the case where the rotation of particles is neglected = , it is straightforward to show that the two expressions become equal.

Second-order work in DEM
The second-order form generated in the REV as a result of implementing the resultant capillary force in DEM reads Similar to the first-order work (Section 4.1), since the current DEM model of unsaturated condition cannot describe the biphasic condition, in Equation (78), we are looking at the second-order work associated with the net stress. In order to obtain variation of , (Equation 4), we first write it in the initial configuration and then calculate its increment: Recalling the relations d = −1 0 d 0 , d = −1 0 d 0 and Equations (63) and (64), Equation (79) can be written in Eulerian form as (virtual) kinematics to the REV and calculating the works of all the internal forces through writing the weak equilibrium equations of different material phases. The obtained micromechanical relations show the existence of additional contributions with respect to the relations presented in the literature. The additional terms in the first-order work expression pertains to the work of fluid-fluid surface tension in stretching the membrane between the fluids and in changing the topology of the intersection line of three phases. In case of the second-order work, the additional contribution comes from the work of suction on changing the topology of the fluid-fluid interface. Future studies should be focused on evaluating the relative importance of the newly proposed terms explained in the previous paragraph using for instance micro-scale numerical simulations. Using such framework, the implications of the new terms on failure prediction of unsaturated granular materials should also be investigated.
It is also important to mention that due to the dependence of the additional terms in the proposed first-and secondorder work expressions on the kinematics of the unsaturated material at the micro-scale, these novel formulations are suitable for use in particle scale models of unsaturated granular media, such as DEM numerical framework or one including coupling between DEM (for solid particles) and Lattice-Boltzmann method for fluids. The benefit of the extended formulations for existing continuum-scale models of unsaturated soils, such as Barcelona Basic Model, 20 is not straightforward, as the state variables in those models do no cover all the features required for our enriched formulations. As such, a re-formulation of those models might be required to introduce new micro-mechanics-based state variables.
Next, we formulated the macroscopic stress of the REV using the first-order work principle. The obtained macroscopic stress coincides with the one derived in Duriez et al. 11 based on volume averaging. This validates the assumptions made in arriving at the current first-and second-order work expressions. An important assumption made was that the interfaces between the three phases do not enter the weak equilibrium equation of the REV as separate material phases.
At the end, some aspects of modeling the unsaturated granular material in DEM were discussed. In particular, the equivalent capillary force, which is currently implemented in DEM for simulating the pendular regime, was shown to produce the same first-and second-order works in the REV as the distributed capillary forces only in the infinitesimal deformation and rotation regime and when the rotation of the particles is neglected. As such, the current DEM simulation of pendular regime can only reproduce the average stress in the REV, as shown by Duriez et al., 11 while it is inadequate in terms of reproducing the internal works of the REV.

A C K N O W L E D G M E N T S
This work was developed within a discovery grant provided by the Natural Sciences and Engineering Research Council of Canada (NSERC). Their funding is gratefully acknowledged.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.