International Journal of Geomechanics

The finite-element formulation of a coupled fluid flow and geomechanics for two-phase fluid flow through fractured porous media are presented. Two porosities, pores and fractures, and five phases are introduced. The two fluids are taken as wetting and nonwetting. The governing equations are derived based on the theory of poroelasticity, the effective stress principle, and the balance equations of mass and momentum, taking into account the solubility of nonwetting fluid into wetting fluid. Spatial and temporal discretization of the governing equations has been realized through the Galerkin method and the finite-difference technique, respectively. A three-dimensional numerical code has been developed and validated based on previously published data. Various applications of the model have been demonstrated through three field-scale examples.

The present study deals with the stress produced in an irregular heterogeneous self-reinforced composite half-space due to a load acting as a shearing force and moving along the free surface. The expression of shear stress produced in this case was obtained in closed form. The effect of heterogeneity, reinforcement, maximum depth of irregularity, and the irregularity factor on shear stress is the highlight of the study. In addition, the effect of different shapes of irregularity on shear stress is discussed through the irregularity factor. Moreover, a comparison is made for three different cases of irregularity, namely, rectangular, parabolic, and no irregularity. As a special case of the problem, the stress produced because of a moving load acting as a shearing force in an isotropic half-space with and without heterogeneity and irregularity is discussed. A comparative study of a self-reinforced composite medium and an isotropic medium in the present context of the problem was performed using graphs, and important peculiarities were traced out.

The knowledge of fluid flow through rock fractures is directly related to hydrocarbon migration, waste disposal, and carbon dioxide sequestration. The hydraulic nature and response of the fractures are directly controlled by the roughness of the fracture surfaces. However, this parameter is hard to understand because it can behave differently under different ambient conditions. The prevalent controlling parameters are the fracture inflow pressure, aperture of the fracture, and shearing displacement during flow. To understand the influence of these parameters, a systematic study was carried out numerically on different fracture geometries. In this paper, two-dimensional fractures with different surface roughness were simulated in a finite-element modeling (FEM) program, and the fluid-flow parameters were evaluated. The Navier–Stokes (NS) equation was used to model the fluid flow through the roughness profiles generated using Barton’s joint roughness coefficient. By simulating the laminar fluid flow through the NS equation and predicting the particle transport using a streamline particle-tracking method, the flow-velocity profiles, outlet-pressure distribution, Reynolds number, shear rates, and particle transmissivity were measured. The parameters at different locations along the length of the fractures were compared to identify changes in the fluid flow. The models show that local undulations have considerable effect on the fluid flow. The velocity and shear-rate evolution, pore-pressure distribution, and Reynolds number of the flow indicate the presence of a strong wall effect on the fluid flow. The aperture and shearing displacements of the fracture walls also have significant control over the wall effect.

In this study, the author considers the propagation of Love waves in an inhomogeneous viscoelastic layer over a gravitational half-space when the upper boundary is assumed to be free. The proposed model is solved to obtain the dispersion relationship for Love waves in the Voigt-type layer under the effect of gravity, inhomogeneity, and initial stress. The relationship between the phase velocity and wave number of Love waves is expressed and then demonstrated with an interactive code. In the absence of various parameters of the medium, the results are in agreement with the classical results. They may be useful for understanding the nature of seismic wave propagation in geophysical applications and in the field of earthquake engineering.

Because most geotechnical analyses may involve elastoplastic geomaterials, a robust solution scheme is of critical importance to the efficiency of the entire finite-element (FE) computation. To accelerate the Krylov subspace iterative methods, some preconditioning techniques have been developed based on the factorization of elastic stiffness. However, the idea of constructing a preconditioner based on elastic stiffness is heuristic. In this study, a class of Neumann expansion preconditioners constructed from the constant (elastic) partition and varying (plastic) partition of the elastoplastic stiffness matrix was proposed. Based on two numerical examples, the performances of truncated Neumann expansion preconditioners were examined with associated and nonassociated soil plasticity considered, respectively. It is interesting to note that the convergence behaviors of truncated Neumann expansion preconditioners closely depended on the approximation to the elastic stiffness part as well as the truncated terms. Numerical experiments also revealed that in most large-size elastoplastic FE analyses, zero-order Neumann expansion preconditioners with appropriate approximations to the elastic stiffness appear to be more appealing, particularly for those applications involving a large amount of nonlinear iterations.

A contact algorithm to determine the aperture evolution of a rock fracture during shearing is reported. The algorithm was based on a theoretical normal closure model and semianalytical dilation model and required only the three-dimensional topography data of the two fracture surfaces at the initial stage before shearing. It allowed the prediction of the aperture distribution of a fracture under normal stress and various shear displacements, which is difficult to observe in coupled shear-flow tests. To ensure the precision of the surface topography data, a laser-scanning profilometer system was used for surface topography measurement. The contact algorithm was then used to predict the aperture evolution of a marble fracture during shearing. The results of aperture variation under normal stress and shear displacement were analyzed and compared with an equivalent aperture derived from back-calculation of flow test results. The results indicated that the contact algorithm is efficient in predicting aperture evolution.

A contact mechanics model, based on the Hertzian elastic contact theory and cementation coating development at particulate scale, was established to predict the time-dependent small-strain stiffness of Class C fly ash–stabilized soils during curing. The cementation coating development model was developed at particulate level based on the Arrhenius law to predict the contact radius growth. A hyperbolic time–temperature relationship was proposed to capture the temperature change of fly ash–stabilized soils and links the pozzolanic reaction rate with curing time. Model-predicted small-strain stiffness was evaluated through both published and experimental test results with good success. The micromechanics modeling indicated that the small-strain stiffness of fly ash–stabilized soil depends on the contact area between fly ash and soil particles and the soil particles’ shear modulus. Most of the small-strain stiffness of the stabilized soil was developed within the first 7 days of curing. In addition, a parametric study and a sensitivity analysis were carried out, which indicated that the proposed contact mechanics model was reliable and robust for predicting the time-dependent small-strain stiffness of soils stabilized with Class C fly ash (or other cementitious stabilizers).

Lightly stabilized granular materials with slow-setting cementitious binders are widely used in the construction of new pavements and rehabilitation of old pavements in Australia and other countries. Higher strength and durability of these materials combined with relatively low cost are reasons they are widely used. The mechanistic-empirical design of pavements with lightly stabilized materials requires key parameters, such as permanent deformation of these materials. A cyclic load-testing facility was developed for studying the permanent deformation characteristics of pavement materials. In particular, a typical pavement structure consisting of a lightly stabilized granular road base material over an expansive soft clay subgrade material was studied using a pavement model testing (PMT) facility. In the PMT facility, the model pavement structure was constructed in a 1.0 × 1.1 × 0.6-m steel tank, and the wheel loading from the vehicle was simulated using a circular steel plate subjected to sinusoidal-type axial loading through an actuator assembly. The results of this study indicate that the permanent deformations of the pavement structure and the subgrade increased with the number of load cycles and were greatly affected by the moisture movement from the subgrade to the stabilized base layer.

Combined vertical and radial consolidation drainage under linear time-dependent loading is presented with varying loading-dependent coefficients of vertical and radial consolidation determined analytically and experimentally. The convolution integral of viscoelasticity theory with variable time was first used to separately predict the average degree of consolidation under time-dependent loading for individual vertical and radial directions. The predicted average degree of combined vertical and radial consolidation was then determined using Carrillo's approach. Four different loading rates were applied in consolidation tests on three types of remolded clay with various plasticity indices to estimate the effects of loading rate on the consolidation settlement. Furthermore, two different spacing ratios of well drain were investigated. The findings indicated that Carrillo's approach for the present three-dimensional consolidation is valid. Furthermore, the viscoelastic settlement prediction in terms of the average degree of consolidation using loading-dependent coefficients of vertical and radial consolidation was more consistent with the experimental results compared with those assuming a constant coefficient of consolidation in the vertical and radial directions. Moreover, specimens with a lower spacing ratio showed a higher average degree of consolidation compared with those with a higher spacing ratio because of the shorter radial drainage path.

A three-dimensional (3D) element partition method for simulating fracture problems is proposed here. In this method, cracks are embedded into elements in a more straightforward manner without any nodal interpolation enrichment or any extra degrees of freedom introduced. When a crack runs across an element, the element is divided into two subbulk elements and one subsurface element. The intersection points between element edges and crack faces are taken as the virtual nodes. Thus, the displacements of the virtual nodes are the extra degrees of freedom. To eliminate these extra degrees of freedom, it is assumed that the displacements of the virtual node are related only to its adjacent nodes at the same side of the crack. The least-square interpolation technique is adopted to characterize their relationship. With this method, the cracked element deformation is related to its neighborhood. The stiffness matrix and nodal force vector are derived. The friction and contact effect between crack faces are implicitly incorporated into the numerical model through the subsurface element. With the 3D element partition method, the crack is allowed to embed into an element without mesh modification or remeshing, which makes the fracture simulation highly efficient.

The prediction of total pull force is critical for static pipe-bursting installation. As a major component of the total pull force, soil expansion prediction is also important; however, there are currently limited methods available for such prediction. In this paper, three cavity expansion solutions were used to predict soil expansion pressure acting upon the expander (bursting head) during static pipe bursting. The determined soil expansion pressure was then used to calculate the expansion force required for pipe bursting with or without consideration of the soil collapse. Calculated results were then compared with those from laboratory pipe-bursting experiments to evaluate the feasibility of the calculation methods. The comparison indicated that the large-strain cavity expansion solution reasonably predicted the soil expansion force. The small-strain solution significantly overestimated the soil expansion force. Additionally, considering soil collapse during pipe bursting resulted in a remarkably higher calculated expansion force than that without consideration; in this study, the soil collapse effect should be taken into account because a sandy soil was used in the laboratory experiments.

The stability on a tunnel face is of great practical significance for safe construction in tunnel engineering. In this work, an advanced three-dimensional (3D) failure mechanism is extended to the face stability analysis of a tunnel driven in anisotropic and nonhomogeneous soils using the kinematic approach of limit analysis. The results obtained from the presented approach are in good agreement with the existing ones, showing that this method can be further applied to the proposed situation in this paper. Two common examples of cohesion variation, a linear variation with depth and a layered soil, are then analyzed. The numerical results show that both anisotropy and nonhomogeneity have a significant impact on the critical face pressure, especially when the cohesion is relatively large or when there is a weaker layer. Therefore, these two factors should be taken into account in tunnel design.

Shallow plate-load tests were simulated using a two-dimensional discrete element procedure together with parameters calibrated from a series of biaxial compression experiments on aluminum rods. The deformation and force transmission characteristics within ground materials were investigated first, and then attentions were focused on studying the influences of the interparticle friction coefficient, the roughness of the loading plate, and the size of loading plate on the testing results, including the ultimate bearing capacity, the deformation behavior, and the failure mode. The well-known scale effect in plate-load tests was also analyzed from the perspectives of both continuum mechanics and discrete element method simulations. It was proven that the nonlinear strength behavior of the materials beneath the foundation is a possible source of the scale effect; however, this effect may be counterbalanced by the boundary restriction effect in laboratory experiments.

The vertical pullout capacity of strip and circular plate anchors embedded horizontally in a layered sandy medium was computed by using the plane strain and axisymmetric lower-bound limit analyses in combination with finite elements and linear optimization. The soil medium below the anchor plate was assumed to be comprised of loose sand. Two different sand layers were considered above the anchor with different combinations of their internal friction angles. For several embedment ratios $(H/B)$, the variations of the pullout factors $F_{\gamma }$ and $F_q$ due to the components of soil unit weight and surcharge, respectively, were computed as a function of $H_{\mathrm{dense}}/H$ for different cases; here, H and $H_{\mathrm{dense}}$ imply (1) the depth of the anchor plate from ground surface and (2) the thickness of the dense sand layer, respectively. The pullout resistance increased continuously with an increase in $H_{\mathrm{dense}}/H$. For a given $H/B$, with the same value of $H_{\mathrm{dense}}/H$, the uplift resistance became greater for a case when the dense sand layer was kept just above the plate rather than placing it close to ground surface. Compared with strip anchors, the uplift factors became considerably greater for circular anchors. The effect of the variation of the unit weights of the two layers above the anchor plate on its pullout resistance was also explored. The results obtained from the analysis compare well with the available theoretical and experimental data.

Knowledge of wave propagation across rock joints is necessary for understanding distinctions between wave attenuation and lithology. Whenever a propagating seismic wave strikes a rock joint at an angle not normal to the joint, wave transformation occurs. The rock joints existing in the field are not oriented in the condition for normal incidence of the wave; hence, wave transformations are very common. In this paper, transformations of P and S_v waves have been modeled numerically with the help of a distinct element method. Validation of the numerical model has been carried out by comparing the wave velocities obtained from bender/extender elements. The dependence of the transformation of obliquely striking waves on the frequency of incident waves and stiffness of the joint has been studied. The velocities of propagating waves and the ratios of stress amplitudes of propagating waves with respect to the incident wave have been obtained and compared for all the possible angles of incidence. The wave velocities and amplitudes of propagating stress waves are found to be dependent on the angle of incidence of the waves at the rock joint. It was observed that the transformed waves that are of the same type as the incident wave contain more energy.

To investigate the vibrational fluidization of granular media, 60 samples of three different granular materials were tested at six normal stresses and a wide range of vibration accelerations. The experiments were conducted on a modified vibrational direct shear apparatus, in which vibration was applied horizontally in the direction of shear. The peak, residual, and vibroresidual shear strength envelopes were obtained for the tested materials at different vibration intensities. The results show that increasing vibration intensity reduces the friction angle of granular materials and increases the normal stress, σ_f, below which granular material is fluidized. The particle shape effect on the strength loss due to vibration was also observed.

The GeneSys project aims to develop a low-permeability, sedimentary reservoir of the Northern German Basin for the provision of geothermal energy to heat the Federal Institute for Geosciences and Natural Resources (BGR) located at the Geo Center in Hanover, Germany. In May 2011, a massive fracturing (frac) operation was performed to create a large heat exchanger at a depth of 3,700 m. Subsequent to the frac operation, two low-rate injection tests were conducted below fracture-closure pressure. The injection phases of the two tests provide evidence for a large and highly conductive fracture, although proppants were not used during the frac operation, indicating that the fracture stays open because of self-propping effects, at least in the explored range of effective stress. Hydraulically active fracture areas of more than 0.5 km² were calculated from an analysis of the pressure transient during the low-rate injection tests relying on matrix permeabilities determined in laboratory experiments. In principle, such areas are sufficient to yield geothermal power in the range of 2 MW (thermal performance). The hydraulically active fracture area and estimated fracture transmissibility diminished between the two low-rate injection tests and, hence, with increasing effective normal stress. The authors suppose this observation is because of a hydromechanical component of self-propping effects. Similarly, the mismatch between the shut-in observations and conventional hydraulic models indicates that hydromechanical processes affect the pressure transient by continuously decreasing the fracture aperture.

A new simplified method is proposed to compute the active earth pressure acting on the backface of a rigid retaining wall undergoing horizontal translation. The effect of soil arching for cohesive backfill soil and friction mobilized along the wall–soil interface was considered. Analytical expressions to determine the slip surface angle and the coefficient of active earth pressure were obtained using the limit-equilibrium approach. These expressions were used for the horizontal flat-element method to obtain theoretical formulas for the active earth pressure, active thrust, and its point of application. Additionally, an implicit solution was derived for the depth of the tension crack in the backfill. A parametric study was undertaken to assess the effects of cohesion, unit weight, friction angle, surcharge pressure, and wall–soil friction angle on the active earth pressure, as well as the effects of the friction angle of backfill soil and wall–soil friction angle on the slip surface angle. The results of the proposed method were then verified against the existing test data and the predictions from the Coulomb and Rankine theories. The results show that the proposed method yielded satisfactory results.

This paper presents a series of unconsolidated-undrained (UU) triaxial compression tests for investigating the behavior and failure mechanism of geotextile-reinforced clay and the effects of sandwiching nonwoven geotextile in a thin layer of sand (sandwich technique) on improving the shear strength of reinforced clay. Test variables include confining pressures, the number of geotextile layers, and thicknesses of the sand layers. The mobilized tensile strain of reinforcements, estimated according to the residual tensile strain by using a digital image–processing technique, was used to directly quantify the effects of soil-geotextile interaction on the shear-strength improvement. The test results showed that the shear strength of reinforced clay increased as the number of geotextile layers was increased. Failure patterns were changed from classical Rankine-type failures for unreinforced soil specimens to bulging (barrel-shaped) failures between adjacent geotextile layers. The effectiveness of reinforcing clay by applying nonwoven geotextile can be attributed to an increase in the apparent cohesion of the reinforced clay specimen. Regarding the sandwich technique, the test results revealed that layers of sand encapsulating the reinforcement can effectively provide an improved soil-geotextile interaction, thereby enhancing the shear behavior of reinforced clay. The shear strength increased as the thickness of the sand layer was increased. An optimal value of sand-layer thickness for maximum shear-strength improvement was not observed at large confining pressures. An appreciable shear-strength improvement was still observed when the sand-layer thickness was increased from 15 to 20 mm at σ₃ > 100 kPa. The sandwich technique contributes to shear-strength improvement by increasing the friction angle of reinforced specimens. The mobilized tensile strain and force of the geotextile increased as the number of geotextile layers, thicknesses of the sand layers, and confining pressure were increased. The mobilized tensile strain and force were strongly correlated to the strength difference between reinforced and unreinforced soil. This experimental finding demonstrated that mobilized tensile strain and force directly contribute to the shear-strength improvement of reinforced clay.

Soil exhibits creep behavior, in which the deformation of soil develops with time at a state of constant effective stress. A series of triaxial compression creep tests was performed on a saturated till to investigate the creep behavior at different effective stresses and strains of the samples. In the tests, the saturated till was allowed to creep in single, multiple, and stepwise stages. The dependency of creep rate of the saturated till on time was examined using the test data. It was found that the time-dependent creep rate models are incapable of describing the creep behavior of soil under the condition of multiple-stage creep, i.e., complex load conditions. A strain-dependent and stress-dependent creep model of the saturated till subject to triaxial compression was proposed. In the model, the creep rate of soil was formulated as a function of inelastic strain and effective stresses instead of time and effective stresses. This model could describe the dependencies of soil creep rate on time, strain, and effective stress observed in the triaxial compression creep tests. Test data obtained in the creep tests were used to verify the proposed model. Moreover, a methodology to accelerate the creep test for soil was developed.

A consolidation analysis method for an earth-rockfill dam based on the coupling of physical state, stress deformation, and seepage was established. In this method, a permeability coefficient model considering physical state and shear stress level was introduced, and it was integrated with Biot’s consolidation theory to simulate the coupling relationship. Taking the Nuozhadu high earth-rockfill dam as an example, a parameter back-analysis method using a neural network and evolutionary algorithm was applied to back-analyze the permeability coefficient model parameters of the core material based on observation data of pore water pressure in the dam. Then, consolidation numerical analyses were performed by using the back-analyzed parameters and the constant permeability measured in situ, respectively. The spatial and temporal distribution of pore water pressure calculated by the back-analyzed parameters coincided quite well with the monitoring data during construction and operation. However, the numerical results with constant permeability are different from the observed values, which indicates that the permeability coefficient model and the consolidation analysis method can be applied effectively, reasonably, and necessarily to numerical simulations for high earth-rockfill dams.

Fine-particle migration and associated internal erosion are major concerns for dam safety. Granular filters are used to prevent fine-particle migration, and several empirical models have been introduced for the design of these filters. Few computational techniques have tracked particle transport through the filters. This paper presents a three-dimensional transient fully coupled pore-scale model used to study the mechanism of fine-particle migration in granular filters. A pore-scale idealization of the fluid was achieved by using the lattice Boltzmann method, and the solid phase was modeled at a microscale using a discrete element method. The fluid forces applied on the particles were calculated on the basis of the momentum exchange between the fluid and particles. The proposed numerical technique was used to model the migration of base-soil particles through granular filters of different particle sizes. Results of conducted simulations provided the erosion percentages and flow rates during the simulations. The general behavior of the filters agreed with observations in laboratory experiments and with filter-design criteria reported in the literature. The proposed computational framework can be used effectively to model fine-particle migration at a pore scale with minimal assumptions, which would add a new dimension to filter-design problems.

Mechanical properties of methane hydrate–bearing soils are influenced by the surrounding temperature and pore pressure. Studies of such influences are of great significance for the safe exploration of methane hydrate. First, a thermo-hydro-mechanical bond contact model is introduced to elucidate the microscopic contact scale behavior of grains with methane hydrate bonds. Second, such a model is incorporated into the distinct-element method (DEM), and a dimensionless temperature-pressure distance parameter is introduced. The influences of temperature and pore pressure along heating and depressurization paths are analyzed by conducting DEM biaxial compression tests, which are also compared with the results of laboratory triaxial compression tests conducted by others. The variation rules of the macromechanical properties (i.e., strength, elastic modulus, peak friction angle, cohesion, and dilatancy angle) with changing temperatures and pore pressures are then analyzed for different effective confining pressures and methane hydrate saturations in the DEM biaxial compression tests, which helps to provide an initial reference for practical applications. Finally, the microscopic mechanisms that affect the macromechanical characteristics are investigated via the ratio of intact bonds under different loading conditions. The results demonstrate that the DEM biaxial compression tests incorporating the bond contact model can efficiently capture the influences of temperature and pore pressure on the macromechanical properties of methane hydrate–bearing soils. It is shown here that temperature and pore pressure affect the mechanical properties of interparticle methane hydrate and, in turn, influence the macroscopic mechanical behavior of methane hydrate–bearing soils.

The aim of this study was to highlight the importance of the numerical approach in realistic support design for tunneling. The Komurhan Tunnel, which is to be constructed between the Elazig–Malatya highway, located east of Elazig City in eastern Turkey, was selected as the application site of this study. The lithology of this area consists of Late Cretaceous–age ophiolitic rocks, known as Komurhan ophiolites. Field and laboratory studies were performed to determine the geotechnical properties of rock mass and intact rock materials. Empirical and numerical approaches were used, and the results were compared, focusing on tunnel design safety. The rock mass was classified in terms of the rock mass rating (RMR), rock mass index (RMi), Austrian tunneling method, and geological strength index rock mass classification systems, and then the support systems for the tunnel were determined accordingly. The support systems obtained from the classification systems were also analyzed using software based on the FEM. On the basis of the numerical analysis, a more reliable and safe support design was recommended. The results show that the support systems suggested by the RMi method were not sufficient for the walls; however, the support systems were found to be safe for the crowns. Furthermore, the support systems suggested by the RMR system did not provide enough stability for both walls and crowns. The support system recommended in this study reduced the thickness of the plastic zone and total displacement values.

The effect of out-of-plane stress on theoretical solutions of stress, displacement, and plastic radius for a circular tunnel were studied by using the generalized Hoek–Brown (HB) failure criterion and nonassociated flow rule in elastic–brittle–plastic rock masses. A technique of using equivalent Mohr–Coulomb and generalized HB strength parameters was adopted to compare results of the developed method and published results. The presented solution was validated by existing theory. The results of the developed method were found to be larger than those previously published, which were determined through programming.

A newly built roller-compacted concrete gravity dam is located on the Yalongjiang River in southwest China. Complex geological conditions greatly puzzled the designers of the project. There are abundant joints within the riverbed that form the double-inclined wedge foundation. This is a typical mode of instability of concrete gravity dams. The equal-K method, which falls into the category of conventional limit-equilibrium methods, is provided in the design codes of gravity dams to evaluate the stability of a dam against sliding along deep failure surfaces. The factor of safety calculated by the equal-K method is determined by the angle of the force acting on the interface of the two inclined wedges, denoted by $\phi$, which is unspecified in the design codes. On the basis of Sarma’s assumption, this study proposes an algorithm for the determination of $\varphi$, which needs only slight improvement of the algorithm in the existing codes. The justification of the proposed procedure is demonstrated by showing that the results from the proposed procedure are very close to those from a rigorous limit-equilibrium method in which all equilibrium conditions are satisfied.

By using the lower-bound finite-element limit-analysis technique, the bearing capacity of a strip footing placed over an embankment was determined with the incorporation of pseudostatic horizontal seismic body forces. The variation of the bearing capacity factors $N_c$, $N_q$, and $N_{\gamma }$ were calculated as the function of the friction angle of soil $\varphi$, slope angle $\beta$, horizontal seismic acceleration coefficient $k_h$, and setback distance $S$. On the basis of the magnitude of $\varphi$, the $S$ values were varied from 0B to 15B. From the results, it can be understood that the bearing capacity factors decrease with the increase in $\beta$ and $k_h$. With the increase in $\varphi$ and $S$, the bearing capacity factors were found to increase. After a particular setback distance, the effect of the embankment slopes was found to be negligible. Failure patterns were also studied.

Darcy’s law of fluid flow and the stress equilibrium equations of anisotropic media, coupled with the mass conservation equations of fluid, are studied on the basis of the model of constitutive equations for both effective stress and porosity of anisotropic saturated soils, in which a new compressible model of soils and anisotropic effects are introduced. This technical note gives a generalized consolidation equation of anisotropic saturated soils. The results show that for isotropic soils, the consolidation equation derived degenerated into Biot’s theory in some approximation, and in the one-dimensional case, there were some subtle differences among Biot’s theory, Terzaghi’s theory, and the current theory. In the final part of this technical note, a set of generalized consolidation equations of orthotropic saturated soils are deduced, and some simple solutions are obtained.