Seepage Force and Its Direct Mechanical Effects in Hydrate- Bearing Porous Media

The direct mechanical effects of seepage force on the behavior of hydrate-bearing porous media (HBPM) have often been neglected in previous studies, which may lead to inaccurate predictions of the mechanical behavior of HBPM under seepage conditions. Here, we propose an extended three-phase physical model for unsaturated HBPM, including gas, capillary water, and generalized solid skeleton (GSS). Based on the model, the force balance equations for the three phases are formulated. Performing a force analysis of generalized solid particles under seepage conditions, we find that the tangential force acting on the generalized solid particle is the seepage force in HBPM. By combining this finding with the formulated balance equations, we derive the expression for the seepage force, which can distinguish the mechanical effects between the tangential force and normal force. The stresses, induced by the tangential force (i.e., seepage force), are divided into two types: one acts on the cross sections of generalized solid particles and the other on the contacts between particles. Neither of them is transmitted through the GSS. The former mainly causes the particles themselves to compress, whereas the latter primarily influences the sliding of the particles at contacts. Based on the mechanisms of these two stresses, their effects on the mechanical behavior of HBPM are quantified, which provides a new insight to evaluate the direct mechanical effects of seepage force.


Introduction
Seepage is an important flow phenomenon in HBPM, and it refers to the flow of water in saturated HBPM or of gas and water in unsaturated HBPM under the gradient of total head [1]. In saturated or unsaturated HBPM, the characteristics of seepage largely determine the magnitude of fluid flow velocity [2], the fluid pressure and saturation distribution [3,4], and the stress and strain distribution of HBPM [5], which further affects the stability of hydrate-bearing sediments [6,7]. As a result, the analysis of the seepage characteristics of HBPM has become an active field in the efficient and safe exploitation of gas hydrate.
Some researchers have performed a large number of experiments by macroscopic and microscopic techniques in the laboratory to investigate the seepage characteristics of HBPM, such as the variation of gas/water relative permeability with different hydrate habits and saturations [8][9][10], the impact of various pressure gradients on gas production rate and water relative permeability [11,12], and the evolution mechanism of relative permeability associated with capillarity and drying-wetting cycles [13,14]. On the other hand, many numerical simulations have been conducted to evaluate the influence of different factors on the gas/water production rate during hydrate exploitation, such as the permeability of host porous media [15,16], the methods used to cause hydrate to dissociate [17][18][19], and the heterogeneity of hydrate saturation [20]. It should be noted that whether the studies were performed by experimental tests or by numerical simulations, they mainly focused on the flow behavior of fluids (i.e., gas and water) in HBPM. The flow behavior, to some extent, reflects the action of the solid phase (i.e., mineral particle and hydrate) in HBPM on the fluids. However, the corresponding reaction is less considered during hydrate exploitation.
The reaction mentioned above represents the force applied by the flowing fluids to the solid phase, usually referred to as seepage force. The mechanical effects of seepage force are generally considered in such a way that the reaction of seepage force results in the variation of pore-fluid pressure, which causes the variation in the effective stress and mechanical response of HBPM [5,[21][22][23]. In fact, the way of considering seepage force only involves its indirect mechanical effects but ignores its direct mechanical effects.
The objective of this work is to identify the mechanism of seepage force in HBPM and to quantify the seepage force and its direct effects on the mechanical behavior of HBPM. For this objective, we first analyze the interactions between hydrates, solid particles, gas, and water and build an extended threephase physical model including the gas phase, capillary water phase, and GSS phase. Based on this model, the independent force balance equations for these three component phases are formulated. Then, a force analysis considering the fluid flow at particle scale is performed, the mechanism of the seepage force is identified, and the expression for the seepage force is derived. Finally, the stresses, induced by the seepage force, acting on the GSS are analyzed, and their direct effects on the mechanical behavior of HBPM are quantified.

An Extended Three-Phase Physical
Model for HBPM HBPM is generally composed of hydrate, gas, capillary water, and solid particles, particularly during the process of hydrate dissociation through depressurization, thermal stimulation, or CO 2 injection [24][25][26][27][28]. The location of hydrate formation depends largely on environmental conditions, such as gas source, water saturation, and the type of host porous media [29,30]. Different environmental conditions lead to different pore-scale habits of hydrate in HBPM. The habits of hydrate are commonly grouped into four categories: cementation, particle coating, pore filling, and load bearing [29,[31][32][33]. The cementing hydrates form at interparticle contacts and bond neighboring particles together. In this case, even a small number of hydrates can remarkably enhance the strength and stiffness of HBPM [34][35][36]; owing to the hydrophilic nature of host porous media, the particle-coating hydrates crystallize on the whole surface of particles. The hydrates of this category can, to a certain extent, bond some particles and thus give rise to a restriction of the movement of these particles [37,38]; the pore-filling hydrates crystallize on the partial surface of particles and enlarge freely into the middle of pores without bonding adjacent particles together. Compared with the cementing hydrates, the pore-filling hydrates are believed to have a relatively small influence on the shear strength and stiffness of HBPM, but they can still restrict the movement of surrounding particles [37,39]; the loadbearing hydrates evolving from the pore-filling hydrates bond particles together into a cluster structure when hydrate saturation is beyond 40% [40]. In particular, the hydrate saturation can be higher than 80% in some local regions where the hydrates may directly bear loads [41,42]. Although hydrates with different pore-scale habits have different contributions to the shear strength and stiffness of HBPM, all of them can be regarded as a structure to bear and transfer loads. Additionally, the interactions between gas and capillary water create gas-water interfaces. The interfaces have the similar mechanical behavior to an elastic contractile skin and, especially those occurring at particle contacts, can pull solid particles together and thus help to increase the strength and stiffness of porous media [43,44]. Hence, the gas-water interfaces can also be regarded as another structure to bear and transfer loads. In addition to the hydrates and gas-water interfaces, the solid particle skeleton in HBPM that is formed by the frictional and interlocking mechanisms among particles is the key structure to bear and transfer loads [45,46]. Considering that the hydrates, gas-water interfaces, and solid particle skeleton all can bear and transfer loads, they could be combined into a generalized structure, called GSS. In such a case, the GSS, gas, and capillary water constitute an extended three-phase physical model ( Figure 1) that can be employed to describe the mechanical behavior of HBPM. For simplicity, the case of particle-coating hydrates is only illustrated. Based on the extended three-phase physical model presented above, a series of definitions can be given to describe the volume-saturation properties of HBPM. For an RVE of HBPM, the following equation must be satisfied ( Figure 2): where V g , V cw , and V gs are the volumes of gas, capillary water, and GSS, respectively, and V is the total volume of RVE. The volume fraction of each phase is defined as the ratio of the volume of the corresponding phase to the total volume of the RVE: where n g , n cw , and n gs are the volume fractions of gas, capillary water, and GSS, respectively, as illustrated in Figure 2.
Combining equations (1), (2), (3), and (4), we have n g + n cw + n gs = 1: The pore spaces occupied by the gas and by the capillary water constitute the effective pore space of HBPM, and the following equation can be given: where V eff v is the effective pore volume. Based on the above definition, the volume fraction of the effective pore space n eff v is then The proportion of the effective pore space occupied by the capillary water is expressed as the effective saturation of the capillary water ( Figure 2): Similarly, the area fraction for each phase is defined as the ratio of the area of the corresponding phase to the total area of the RVE. In general, the area fraction is postulated to be equal to the volume fraction in the case of homogeneous porous media [47].
Based on the analysis of the pore-scale habits of hydrate and the interactions between the gas and capillary water, an extended three-phase physical model for HBPM is well established, as shown in Figure 1. The model provides a physical basis for the derivation of the force balance equations for three individual phases in HBPM.

Force Balance Equations for Three Individual Phases in HBPM
The force balance equations for three individual phases in HBPM are the prerequisites for determining the seepage force caused by the flowing fluids and its direct effects on the mechanical behavior of HBPM. The force balance equations for individual phases in unsaturated soils have been formulated by the differential element-based force balance method [48][49][50]. This method is based on the widely accepted concepts of force balance and can explicitly take the interaction forces between different phases into account. Therefore, this method will be adopted to formulate the force balance equations for the gas, capillary water, and GSS contained in HBPM.

Force Balance Equation for the Gas Phase.
Formulating the force balance equation for the gas phase needs to calculate the surface forces, gravitational force, and interaction force that act on the differential element of the gas phase. The surface forces are calculated as the gas pressure times the area occupied by the gas. The gravitational force is equal to the product of the unit weight of the gas and the corresponding volume. The interaction force is expressed as the interaction force per unit volume times the corresponding volume. For simplicity, only the forces in the y-direction are shown in Figure 3. The equilibrium condition for the differential element of the gas phase requires that the resultant force should vanish. Summing the force components in the y-direction yields the balance equation for the gas phase: where u g is the gas pressure, f gs gy is the interaction force per unit volume between the gas and the GSS; ρ g is the density of gas, and g is the gravitational acceleration.

Force Balance Equation for the Capillary Water Phase.
In analogy with the force analysis of the gas phase, the force components acting on the differential element of the capillary water phase in the y-direction are shown in Figure 4. Capillary water Generalized solid skeleton Figure 2: Schematic illustration of the volume-saturation relations of an RVE of HBPM.

Geofluids
The summation of the force components in the y-direction gives the balance equation for the capillary water phase: where u w is the capillary water pressure, f gs cwy is the interaction force per unit volume between the capillary water and the GSS, and ρ w is the density of water.

Force Balance Equation for
the GSS Phase. Compared with the differential elements of the gas phase and capillary water phase that are only subjected to one type of surface force, the differential element of the GSS phase is subjected to two types of surface force. One represents the forces exerted by the adjacent differential elements of the GSS phase, and the other refers to the forces induced by the gas pressure and capillary water pressure. As the gas and capillary water together occupy the effective pore space, a difficult problem encountered is how to determine the magnitudes and corresponding areas of the stresses acting on the differential element of the GSS induced by the second type of surface force.
To solve the above problem, we shall first consider a saturated system composed of two idealized spherical solid particles coated with hydrates in equilibrium under hydrostatic pressure, as illustrated in Figure 5.
According to the equilibrium condition for section 1-1 ( Figure 5(b)), we can obtain where σ w sp is the stress acting on the cross section of the solid particle induced by the hydrostatic pressure and A sp is the cross-sectional area of the solid particle.
Similarly, using the equilibrium condition for section 2-2 ( Figure 5(c)), we can get where σ w c is the stress acting on the contact between particles induced by the hydrostatic pressure and A c is the contact area.
In analogy with the saturated system, an unsaturated system is analyzed that includes two idealized spherical generalized solid particles (composed of solid particles, particlecoating hydrates, and gas-water interfaces) in equilibrium under the gas pressure and capillary water pressure, as shown in Figure 6. In order to apply the approach developed in the saturated system to the unsaturated system, a basic assumption analogous to that adopted by the theory of mixtures [51] is made here. It is assumed that the gas and capillary water independently fill the total pore space of the unsaturated system according to their respective volume fractions. In this case, this system can be regarded as two subsystems. One refers to a system where the pore space is fully filled with the gas under a homogenized gas pressure of ½n g /ðn cw + n g Þu g (Figure 6(b)); the other corresponds to a system where the capillary water under a homogenized water pressure of ½n cw /ðn cw + n g Þu w completely occupies the pore space ( Figure 6(c)). From Figure 6(b), on the basis of the equilibrium condition of section 3-3, the following equation can be obtained: where σ g sp denotes the stress acting on the cross section of the generalized solid particle induced by the homogenized gas pressure.

Geofluids
Likewise, according to the equilibrium condition of section 4-4, the following equation can be obtained: where σ g c denotes the stress acting on the contact between generalized solid particles induced by the homogenized gas pressure.
From Figure 6(c), using the equilibrium conditions of sections 5-5 and 6-6, equations (15) and (16) can be given, respectively: where σ cw sp denotes the stress acting on the cross section of the generalized solid particle induced by the homogenized capillary water pressure.
where σ cw c denotes the stress acting on the contact between generalized solid particles induced by the homogenized capillary water pressure.
Extending the analysis results obtained in the unsaturated system to a differential element of the GSS in HBPM, we can easily compute the magnitudes and corresponding areas of the stresses acting on the differential element induced by the homogenized gas pressure and by the homogenized capillary water pressure. These stresses are illustrated in Figure 7. Besides, Figure 7 also illustrates the gravitational force, interaction forces per unit volume, and the stresses exerted by adjacent differential elements.
Summing the force components in the y-direction yields the balance equation for the GSS phase: where τ xy ′ is the shear stress acting on the x-plane in the y -direction, σ y ′ is the normal stress acting on the y-plane, τ zy ′ is the shear stress acting on the z-plane in the y-direction, f g gsy is the interaction force per unit volume between the GSS and the gas, f cw gsy is the interaction force per unit volume between the GSS and the capillary water, and ρ gs is the density of the solid particle. The detailed derivation of equation (17) is presented in Appendix A.

Identification and Quantification of Seepage Force in HBPM
When the gas and capillary water flow in the pore space of HBPM, they are subjected to the seepage resistance exerted by the GSS. On the contrary, the GSS are subjected to the seepage force imposed by the flowing gas and capillary water. To identify the seepage resistance and seepage force in unsaturated HBPM, it is necessary to make an identification of the seepage resistance and seepage force in saturated HBPM.
4.1. Seepage Force in Saturated HBPM. A saturated system, which is composed of two idealized spherical solid particles coated with hydrates, is used to identify the seepage resistance and seepage forces under saturated seepage conditions, as illustrated in Figure 8. It can be seen from Figure 8(a) that the solid particles are subjected to the normal and tangential forces exerted by the flowing water. Conversely, the flowing water is also subjected to the normal and tangential forces exerted by the solid particles (Figure 8(b)). The normal forces acting on the flowing water generate pore-water pressure and thus increase the pressure potential of pore water, whereas the tangential forces acting on the water hinder the flow of water and thus reduce the pressure potential of pore water. As such, the tangential forces acting on the flowing water are the seepage resistance exerted by the solid particles. On the other hand, the normal forces acting on the solid particles cause the volumetric strain of particles, but the tangential forces acting on the solid particles promote the movement of particles. Hence, the tangential forces acting on the solid particles are the seepage forces exerted by the flowing water. Since the tangential forces acting on the flowing water and on the solid particles have the same value but opposite directions, the seepage resistance and the seepage force are numerically equal but opposite in direction.
Extending the results from a particle-level analysis to a differential element-level analysis allows us to calculate the seepage resistance and seepage force in saturated HBPM. In this case, according to equation (10), the force balance equation for the water phase in saturated HBPM can be given by where f gs wy is the interaction force per unit volume between the water and the GSS.
When the water flows in the pore space of the differential element of saturated HBPM, it is subjected to the normal and tangential forces exerted by the GSS. The effects of the normal and tangential forces acting on the water phase correspond to the first and second terms on the left side of equation (18), respectively. Therefore, the interaction force per unit volume f gs wy acting on the water phase is the seepage resistance in saturated HBPM. By substituting the relation n cw u w (n cw + n g ) n g u g (n cw + n g ) n cw u w (n cw + n g ) Figure 6: An unsaturated system composed of two idealized spherical generalized solid particles under static gas and capillary water pressure, including two subsystems: (a) an unsaturated system where the effective pore space is occupied by the gas and capillary water together; (b) a subsystem where the effective pore space is fully filled with the gas under homogenized gas pressure; (c) a subsystem where the effective pore space is completely occupied by the capillary water under homogenized capillary water pressure. 6 Geofluids between pore-water pressure and pressure head, u w = ρ w gh p , equation (18) can be rewritten as where h p is the pressure head of water and y is the elevation head. Then, the seepage resistance in saturated HBPM in the y -direction j gs wy can be expressed as where i y is the gradient of the total head of water in the y -direction. Similarly, based on equation (17), the force balance equation for the GSS phase in saturated HBPM can be obtained: where f w gsy is the interaction force per unit volume between the GSS and the water and n gs = 1 − n eff v .When the water flows, the GSS is subjected to the normal and tangential forces exerted by the flowing water. The effects of the normal and tangential forces acting on the GSS associate with the fourth and fifth terms on the left side of equation (21), respectively. Hence, the interaction force per unit volume f w gsy acting on the GSS is the seepage force in saturated HBPM. As demonstrated earlier, the seepage force and the seepage resistance are numerically equal, and the seepage force in saturated HBPM can therefore be expressed as where j w gsy is the seepage force in the y-direction acting on the GSS.
Generalizing the result obtained in the y-direction to the xand z-directions, we have where j w gs is the seepage force vector acting on the GSS exerted by the flowing water and i is the gradient of the total head of water.  7 Geofluids unsaturated HBPM, we can compute the seepage resistance and seepage force in unsaturated HBPM. According to equations (9) and (10), the seepage resistance in unsaturated HBPM can be given by equations (24) and (25), respectively: j gs gy = n g ρ g g ∂ h p-g + y À Á ∂y = n g ρ g gi gy , ð24Þ

Seepage Force in
where j gs gy is the seepage resistance in the y-direction acting on the gas phase, h p-g is the pressure head of gas and h p-g = u a /ρ a g, and i gy is the gradient of the total head of gas in the y-direction. j gs cwy = n cw ρ w g where j gs cwy is the seepage resistance in the y-direction acting on the capillary water phase, h p-cw is the pressure head of capillary water and h p-cw = u w /ρ w g, and i cwy is the gradient of the total head of capillary water in the y -direction.
When both the gas and the capillary water flow, the GSS is subjected to the normal and tangential forces exerted by the flowing gas and capillary water. The effects of the normal and tangential forces acting on the GSS exerted by the gas correspond to the fourth and sixth terms on the left side of equation (17), respectively. In equation (17), the fifth and seventh terms on the left side are attributed to the effects of the normal and tangential forces acting on the GSS exerted by the capillary water, respectively. Therefore, the interaction forces per unit volume, f g gsy and f cw gsy , acting on the GSS are the seepage forces exerted by the flowing gas and by the flowing capillary water, respectively. Performing similar analyses in the xand z-directions and then using the conclusion drawn above that the seepage force is numerically equal to the seepage resistance, we can give the seepage forces in unsaturated HBPM as follows: where j g gs is the seepage force vector acting on the GSS exerted by the flowing gas and i g is the gradient of the total head of gas.
where j cw gs is the seepage force vector acting on the GSS exerted by the flowing capillary water and i cw is the gradient of the total head of capillary water.

Mechanical Effects of the Seepage Force in HBPM
The mechanical effects of the seepage force in HBPM depend on the transfer mechanism of the stresses acting on the GSS induced by the seepage force and the magnitude of these stresses. In order to identify the transfer mechanism, we need to turn to the effective stress equation for HBPM.

Effective Stress Equation for
Unsaturated HBPM. The effective stress equation for unsaturated HBPM can be derived based on four force balance equations for the gas phase, the capillary water phase, the GSS phase, and the whole HBPM. As the first three force balance equations have been formulated, i.e., equations (9), (10), and (17), the total force balance equation for the whole HBPM will be formulated next. Taking a differential element from HBPM and then analyzing the forces acting on it in the y-direction (Figure 9), we can obtain the total balance equation for the whole HBPM: where τ xy is the total shear stress acting on the x-plane in the y-direction, σ y is the total normal stress acting on the y-plane, τ zy is the total shear stress acting on the z-plane in the y-direction, and ρ is the density of HBPM and ρ = n g ρ g + n cw ρ w + n gs ρ gs . The superposition of equations (9), (10), and (17) for the gas phase, capillary water phase, and GSS phase, ∂τ zy ′ ∂z + n g ρ g g + n cw ρ w g + n gs ρ gs g = 0: ð29Þ Considering that both the gas and the capillary water are unable to carry shear stress, the identities τ xy ′ = τ xy and τ zy ′ = τ zy hold. In this case, equation (29) can be rewritten as ∂τ xy ∂x + ∂σ y ′ ∂y + ∂ ∂y n g n cw + n g u g ! + ∂ ∂y n cw n cw + n g u w ! + ∂τ zy ∂z + n g ρ g g + n cw ρ w g + n gs ρ gs g = 0: ð30Þ Comparing equation (30) with equation (28) and then generalizing the result obtained in the y-direction to the xand z-directions, we get σ = σ ′ + n g n cw + n g u g + n cw n cw + n g u w : By substituting equations (2), (3), (6), and (8) into equation (31), the effective stress equation for unsaturated HBPM can be derived: where σ ′ is the effective stress, σ is the total stress, ½S eff cw u w + ð1 − S eff cw Þu g is the neutral stress, and S eff cw is the effective saturation of the capillary water, which varies from 0 to 1.
Many studies [49,50,[52][53][54] have shown that the effective stress is transmitted through the solid skeleton and largely controls the shear strength and volumetric strain behavior of soils. This statement is conducive to the identification of the transfer mechanism of the seepage force.

Mechanical Effects of the Seepage Force in Saturated
HBPM. According to the balance equation (equation (21)) for the GSS in saturated HBPM, it is found that there are three stresses acting on the GSS, i.e., the effective stress, the stress induced by normal forces, and the stress due to tangential forces (or seepage forces). As the effective stress is transmitted through the GSS and governs largely the mechanical behavior of HBPM, the stresses induced by normal forces and by seepage forces are not transmitted through the GSS. After identifying the transfer mechanism of the stress induced by seepage forces, the magnitude of this stress needs to be determined. We now consider a saturated system consisting of two idealized spherical solid particles coated with hydrates under saturated seepage conditions, as illustrated in Figure 10. In this figure, only the seepage forces are illustrated for simplicity because only the mechanical effects of the seepage forces will be analyzed.
For a particle-scale system (Figure 10), the seepage forces acting on the solid particle can be combined into a resultant force. By considering equations (18) and (23), this resultant force can be computed as where F w-s gs is the resultant force acting on the solid particle induced by the seepage forces and A is the total area of the system (Figure 10(a)). In this system, n eff v = ðA − A sp Þ/ A = 1 − ðA sp /AÞ, and A sp is the cross-sectional area of the solid particle.
According to the equilibrium condition for section 1-1 (Figure 10(b)), we get where σ w-s sp is the stress acting on the cross section of the solid particle induced by the seepage forces.
As σ w-s sp is not transmitted through the GSS, it only produces the volumetric strain of the solid particle. Then,  Figure 9: Force components acting on the differential element of the whole HBPM in the y-direction. 9 Geofluids the contribution of σ w-s sp to the volumetric strain of saturated HBPM can be calculated as where ε w-s v-sp is the volumetric strain contributed by σ w-s sp induced by the seepage forces under saturated seepage conditions and C sp is the compressibility of the solid particle. Likewise, using the equilibrium condition for section 2-2 ( Figure 10(c)), we have where σ w-s c is the stress acting on the contact between solid particles induced by the seepage forces and A c is the contact area.
Owing to σ w-s c not being transmitted through the GSS, it only affects the sliding of solid particles. Then, the contribution of σ w-s c to the shear strength of saturated HBPM can be computed as where τ w-s f -c is the shear strength contributed by σ w-s c induced by the seepage forces under saturated seepage conditions, a c is the contact area ratio and a c = A c /A, and ψ i is the intrinsic friction angle of the solid particle.

Mechanical Effects of the Seepage Force in Unsaturated
HBPM. The mechanical effect of seepage force in unsaturated HBPM can be quantified by employing the approach developed in saturated HBPM. We now consider an unsaturated system consisting of two idealized spherical generalized solid particles under unsaturated seepage conditions, as illustrated in Figure 11. In order to employ the approach developed in saturated HBPM, an assumption analogous to that made in Figure 6 needs to be introduced here. In this case, the unsaturated system can be considered two subsystems. One represents a system saturated with gas under a homogenized gas pressure of ½n g /ðn cw + n g Þu g (Figure 11(b)); the other refers to a system saturated with capillary water under a homogenized water pressure of ½n cw /ðn cw + n g Þu w (Figure 11(c)).
In the light of the results obtained in the saturated system, we can easily compute the seepage forces acting on the generalized solid particle in the unsaturated system. For a subsystem fully saturated with gas ( Figure 11(b)), the seepage forces acting on the particle can be composed into a resultant force. This resultant force can be expressed as where F g-s gs is the resultant force acting on the particle induced by the seepage forces in a subsystem fully saturated with gas.
From Figure 11(b), according to the equilibrium conditions of sections 3-3 and 4-4, the following two equations can be obtained: where σ g-s sp and σ g-s c are the stresses acting on the cross section of the generalized solid particle and on the contact between particles, respectively, induced by the seepage forces in a subsystem fully saturated with gas.
Similarly, for a subsystem completely saturated with capillary water (see Figure 11(c)), the seepage forces acting on the generalized solid particle can be combined into a resultant force. This resultant force can be expressed as where F cw-s gs is the resultant force acting on the particle induced by the seepage forces in a subsystem completely saturated with capillary water.
From Figure 11(c), using the equilibrium conditions of sections 5-5 and 6-6, equations (41) and (42) can be given, respectively: σ cw-s where σ cw-s sp and σ cw-s c denote the stresses acting on the cross section of the generalized solid particle and on the contact between particles, respectively, induced by the seepage forces in a subsystem completely saturated with capillary water.
As σ g-s sp and σ cw-s sp are not transmitted through the GSS, they only produce the volumetric strain of the solid particle. Then, the contribution of σ g-s sp and σ cw-s sp induced by the seepage forces to the volumetric strain of unsaturated HBPM can be calculated as n g u g (n cw + n g ) n cw u w (n cw + n g ) n g u g (n cw + n g ) n cw u w (n cw + n g ) Figure 11: An unsaturated system composed of two idealized spherical generalized solid particles under unsaturated seepage conditions, including two subsystems: (a) an unsaturated system where the effective pore space is occupied by the gas and capillary water together; (b) a subsystem fully saturated with gas under homogenized gas pressure; (c) a subsystem completely saturated with capillary water under homogenized capillary water pressure.

Geofluids
where ε us-s v-sp is the volumetric strain contributed by ðσ g-s sp + σ cw-s sp Þ that is induced by the seepage forces under unsaturated seepage conditions.
Because of σ g-s c and σ cw-s c not being transmitted through the GSS, they only influence the sliding of generalized solid particles. Then, the contribution of σ g-s c and σ cw-s c induced by the seepage forces to the shear strength of unsaturated HBPM can be computed as where τ us-s f -c is the shear strength contributed by ðσ g-s c + σ cw-s c Þ that is induced by the seepage forces under unsaturated seepage conditions.

Conclusions
By analyzing the interactions between hydrates, solid particles, gas, and water in HBPM, an extended three-phase physical model is presented, including the gas phase, capillary water phase, and GSS phase. Based on the model, the independent force balance equations for these three phases are formulated, which provide a sound theoretical basis for the derivation of seepage force.
Under saturated or unsaturated seepage conditions, the forces acting on the fluids and on the solid particles are composed of normal and tangential forces. For the fluids, the normal forces generate pore-fluid pressure and thus increase the pressure potential of fluids, whereas the tangential forces restrict the flow of fluids and thus reduce the pressure potential of fluids. Hence, the tangential forces acting on the fluids are the seepage resistance. On the contrary, for the generalized solid particles, the normal forces produce the volumetric strain of particles, while the tangential forces drive the particles to move. Therefore, the tangential forces that act on the generalized solid particles are the seepage forces. According to this finding, for saturated HBPM, the seepage force is expressed as equation (23); for unsaturated HBPM, the seepage forces induced by the gas and capillary water are calculated by equations (26) and (27), respectively. Compared with the equations for seepage force derived by previous researchers [55,56], the equations obtained in this study distinguish the mechanical influence of the tangential force from that of the normal force, which ensures that the seepage force and seepage resistance are equal in magnitude but opposite in direction.
The seepage forces cause the stresses acting on the cross section of the generalized solid particle and on the contact between particles. Based on the effective stress equation (equation (32)) and the balance equations (equations (21) and (17)) for the GSS, it is found that the stresses acting on the cross section of the generalized solid particle and on the contact between particles are not transmitted through the GSS. The former mainly produces the volumetric strain of particles, whereas the latter primarily influences the sliding of the particle at contacts. Therefore, the mechanical effects of the stresses caused by the seepage forces can be evaluated by equations (35) and (37) for saturated HBPM and by equations (43) and (44) for unsaturated HBPM. These equations allow us to make a direct evaluation of the mechanical effects of seepage force.

Appendix Force Balance Equation for the GSS Phase
The equilibrium condition for the differential element of the GSS phase demands that the resultant force should vanish, as shown in Figure 7. Summing the force components in the y -direction yields − σ y ′ + ∂σ y ′ ∂y dy ! dxdz + σ y ′dxdz − τ xy ′ + ∂τ xy ′ ∂x dx ! dydz + τ xy ′ dydz − τ zy ′ + ∂τ zy ′ ∂z dz ! dxdy + τ zy ′ dxdy − σ g sp-U a g sp dxdz + σ g sp-L a g σπ dxdz − σ g c-U a g c dxdz + σ g c-L a g c dxdz − σ cw sp-U a cw sp dxdz + σ cw sp-L a cw sp dxdz − σ cw c-U a cw c dxdz + σ cw c-L a cw c dxdz − f g gsy dxdydz − f cw gsy dxdydz − n gs ρ gs gdxdydz = 0, ðA:1Þ where σ g sp-U and σ g c-U denote the stresses acting on the cross sections of generalized solid particles and on the contacts between generalized solid particles of the upper surface of the differential element, respectively, induced by the homogenized gas pressure (see Figure 7), σ g sp-L and σ g c-L denote the stresses acting on the cross sections of generalized solid particles and on the contacts between generalized solid particles of the lower surface, respectively, induced by the homogenized gas pressure, a g sp is the corresponding area of σ g sp-U and σ g sp-L , a g c is the corresponding area of σ g c-U and σ g c-L , σ cw sp-U and σ cw c-U denote the stresses acting on the cross sections of generalized solid particles and on the contacts between generalized solid particles of the upper surface of the differential element, respectively, induced by the homogenized capillary water pressure, σ cw sp-L and σ cw c-L denote the stresses acting on the cross sections of generalized solid particles and on the contacts between generalized solid particles of the lower surface, respectively, induced by the homogenized capillary water pressure, a cw sp is the corresponding area of σ cw sp-U and σ cw sp-L , and a cw c is the corresponding area of σ cw c-U and σ cw c-L .