Study on Two-Phase Fluid-Solid Coupling Characteristics in Saturated Zone of Subgrade Considering the E ﬀ ects of Fine Particles Migration

: The ﬂuid seepage in saturated zone of subgrade promotes the migration of ﬁne particles in the ﬁller, resulting in the change of pore structure and morphology of the ﬁller and the deformation of solid skeleton, which a ﬀ ects the ﬂuid seepage characteristics. Repeatedly, the muddy interlayer, mud pumping, and other diseases are ﬁnally formed. Based on the theory of two-phase seepage, the theory of porous media seepage, and the principle of e ﬀ ective stress in porous media, a two-phase ﬂuid-solid coupling mathematical model in saturated zone of subgrade considering the e ﬀ ects of ﬁne particles migration is established. The mathematical model is numerically calculated with the software COMSOL Multiphysics ® . The two-phase seepage characteristics and the deformation characteristics of the solid skeleton in saturated zone of the subgrade are studied. The research results show that the volume fraction of ﬁne particles ﬁrst increases then decreases and ﬁnally becomes stable with the increase of time, due to the continuous erosion and migration of ﬁne particles in saturated zone of the subgrade. The volume fraction of ﬁne particles for the upper part of the subgrade is larger than that for the lower part of the subgrade. The porosity, the velocity of ﬂuid, the velocity of ﬁne particles, and the permeability show a trend of increasing ﬁrst and then stabilizing with time; the pore water pressure has no signiﬁcant changes with time. The vertical displacement increases ﬁrst and then decreases slightly with the increase of time, and ﬁnally tends to be stable. For the ﬁller with a larger initial volume fraction of ﬁne particles, the maximum value of the volume fraction of ﬁne particles caused by ﬂuid seepage is larger, and the time required to reach the maximum value is shorter. It can be concluded that the volume fraction of ﬁne particles in the subgrade ﬁller should be minimized on the premise that the ﬁller gradation meets the requirements of the speciﬁcation in actual engineering. Author Contributions: Conceptualization, Y.D. and Y.J.; methodology, J.-s.Z., X.-b.C., and X.W.; software, F.M.; validation, Y.D. and Y.J.; formal analysis, Y.D.; investigation, Y.J.; resources, J.-s.Z.; data curation, Y.J.; writing—original draft preparation, Y.D.; writing—review and editing, Y.D.; visualization, Y.D.; supervision and project administration, J.-s.Z.; funding acquisition, read and


Introduction
Filler for railway subgrade is mostly coarse-grained soil materials, mixed with a certain proportion of fine-grained soil, and has complex pore structure characteristics. When atmospheric precipitation, surface water, or snow melt infiltrate into subgrade and fail to be discharged in time, the local area of the subgrade will be saturated. Under the action of external loads, excess pore water pressure will be generated in saturated zone of the subgrade, and the dissipation of excess pore water pressure mixture. On the basis of fluid mechanics and soil mechanics, W. Chen [32] derived the equilibrium equation of porous media and the continuity equation of pore fluid considering the coupling of stress field and seepage field, and established a fluid-solid coupling model. However, the above researches on the characteristics of fine particles migration and fluid-solid coupling are mostly separate, and seldom involve the influence of fine particles migration on the hydraulic characteristics of coarse-grained soils. There is a lack of systematic research on the coupling effect between fine particles migration, fluid seepage, and deformation characteristics of solid skeleton.
Therefore, based on the theory of two-phase seepage, the theory of porous media seepage, and the principle of effective stress in porous media, this paper establishes a two-phase fluid-solid coupling mathematical model considering the effects of fine particles migration. The mathematical model is numeralized by the software COMSOL Multiphysics ® . The characteristics of fine particles migration, the two-phase seepage and the deformation of solid skeleton are studied.

Basic Assumption
Filler for railway subgrade is mostly coarse-grained soil materials, mixed with a certain proportion of fine-grained soil, and has complex pore structure characteristics. Therefore, the saturated coarse-grained filler can be regarded as a saturated porous medium composed of solid medium (coarse particles and fine particles), fluid medium in pores, and fine particles that erode into the fluid medium that can migrate with the fluid (hereinafter referred to as fine particles). The structure of saturated porous media is shown in Figure 1. In order to establish a two-phase fluid-solid coupling mathematical model in saturated zone of subgrade considering the effects of fine particles migration, the following assumptions are made: (1) The coarse-grained filler is uniform and isotropic linear elastomer, and satisfies the assumption of small strain; (2) The seepage of water in the pores is caused by the pore water pressure gradient and follows Darcy's law; (3) The influence of temperature and chemical action on the properties of subgrade filler is not considered; (4) Solid particles and water are incompressible; (5) The movable fine particles are spheres with the same radius (d p ), and the influence of fine particles on the fluid seepage characteristics and the force between fine particles are ignored.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 24 mechanics, W. Chen [32] derived the equilibrium equation of porous media and the continuity equation of pore fluid considering the coupling of stress field and seepage field, and established a fluid-solid coupling model. However, the above researches on the characteristics of fine particles migration and fluid-solid coupling are mostly separate, and seldom involve the influence of fine particles migration on the hydraulic characteristics of coarse-grained soils. There is a lack of systematic research on the coupling effect between fine particles migration, fluid seepage, and deformation characteristics of solid skeleton.
Therefore, based on the theory of two-phase seepage, the theory of porous media seepage, and the principle of effective stress in porous media, this paper establishes a two-phase fluid-solid coupling mathematical model considering the effects of fine particles migration. The mathematical model is numeralized by the software COMSOL Multiphysics ® . The characteristics of fine particles migration, the two-phase seepage and the deformation of solid skeleton are studied.

Basic Assumption
Filler for railway subgrade is mostly coarse-grained soil materials, mixed with a certain proportion of fine-grained soil, and has complex pore structure characteristics. Therefore, the saturated coarse-grained filler can be regarded as a saturated porous medium composed of solid medium (coarse particles and fine particles), fluid medium in pores, and fine particles that erode into the fluid medium that can migrate with the fluid (hereinafter referred to as fine particles). The structure of saturated porous media is shown in Figure 1. In order to establish a two-phase fluid-solid coupling mathematical model in saturated zone of subgrade considering the effects of fine particles migration, the following assumptions are made: (1) The coarse-grained filler is uniform and isotropic linear elastomer, and satisfies the assumption of small strain; (2) The seepage of water in the pores is caused by the pore water pressure gradient and follows Darcy's law; (3) The influence of temperature and chemical action on the properties of subgrade filler is not considered; (4) Solid particles and water are incompressible; (5) The movable fine particles are spheres with the same radius ( p d ), and the influence of fine particles on the fluid seepage characteristics and the force between fine particles are ignored.   (1) Basic definitions As shown in Figure 1, the micro-model for coarse-grained soil is composed of three parts, which includes solid medium ( s, dm s is mass of the solid and dV s is volume of the solid), fluid medium (l, dm l is mass of the solid and dV l is volume of the solid), and fine particles (p, dm p is mass of the solid and dV p is volume of the solid).
Porosity of porous media is defined as: where: φ is the porosity of the micro-model; dV V is the pore volume of the micro-model (m 3 ); dV is the total volume of the micro-model (m 3 ). The volume fraction of fine particles is: where: C is the volume fraction of fine particles.
(2) Equations of mass conservation According to the law of mass conservation of fluid, when considering the two-phase (fluid and fine particles) seepage in porous media, the mass conservation equations of fluid phase and fine particles phase can be expressed as: where: ρ l is the inherent density of the fluid (kg/m 3 ); ρ p is the inherent density of fine particles (kg/m 3 ); v ls is the velocity vector of the fluid relative to the solid (m/s); q l is the source (sink) item of the fluid (kg/m 3 s); v ps is the velocity vector of fine particles relative to solid (m/s); . m is the mass rate of fine particles migration (kg/m 3 s); ∇ = ∂ ∂x i + ∂ ∂y j + ∂ ∂z k is Hamiltonian operator. (3) Equations of motion The seepage of fluid conforms to Darcy's law, as follows: where: v l is the velocity vector of the fluid (m/s); v s is the velocity vector of solid (m/s); k is the permeability of porous media (m 2 ); p l is the pore pressure of the fluid (Pa); µ l is the dynamic viscosity coefficient of the fluid (Pa · s); g is the gravitational acceleration vector (N/kg). The force of fine particles in liquid is complex. The fine particles is mainly affected by effective gravity, drag force, additional mass force, pressure gradient force, Basset force, Magnus force, Saffman force, and other force, among which additional mass force and pressure gradient force, Basset force, Magnus force, Saffman force have relatively small influence on the movement of fine particles and can be ignored. Therefore, considering the effect of effective gravity and drag force on fine particles, the schematic diagram of the force on a single fine particle is shown in Figure 2.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 24 particles and can be ignored. Therefore, considering the effect of effective gravity and drag force on fine particles, the schematic diagram of the force on a single fine particle is shown in Figure 2. We use r v to represent the difference between the velocity vector of the fluid and the velocity vector of the fine particles ( p v ). Assuming that the velocity of the fine particles differs little from that of fluid in the horizontal direction, and the velocity difference in the vertical direction can be expressed as: When a single particle is in force equilibrium in the fluid, there are: then:

4( )
where, p d is the diameter of fine particles ( m); D C is the coefficient of drag force, which is related to the Reynolds number of fine particles.
According to Equations (6)~(10), the motion equations of fine particles can be obtained: The schematic diagram of the force on a single fine particle.
We use v r to represent the difference between the velocity vector of the fluid and the velocity vector of the fine particles (v p ). Assuming that the velocity of the fine particles differs little from that of fluid in the horizontal direction, and the velocity difference in the vertical direction can be expressed as: When a single particle is in force equilibrium in the fluid, there are: then: where, d p is the diameter of fine particles (m);C D is the coefficient of drag force, which is related to the Reynolds number of fine particles.

Equations of Stress Field
According to the principle of effective stress in saturated porous medium, there are: where, σ ij is the total stress tensor; σ ij is the effective stress tensor; δ ij is the Kronecker symbol, Appl. Sci. 2020, 10, 7539 6 of 22 According to the assumption that the solid skeleton is uniform and isotropic linear elastomer, and satisfies the assumption of small strain, its constitutive equation is: (13) where, ε V is the volumetric strain of the solid skeleton, ε V = ε x + ε y + ε z , and ε x , ε y , ε z is the axial strain in the three directions of x, y, and z, respectively; ε ij is the strain of the solid skeleton; λ, µ is the constant of Lamé, where, E is Young's elastic modulus; v is Poisson's ratio. The displacement of the solid skeleton is expressed in terms of u s , and the relation between strain and displacement is: Then, the stress balance equation of porous media can be expressed as: where, R i is the volume force of the porous medium (N); where, ρ is the average density of the porous medium (kg/m 3 ); ρ s is the inherent density of the solid (kg/m 3 ). Then, the governing equation of the stress field of porous media can be obtained from Equations (12), (13), and (16)~ (18).

Migration Equation of Fine Particles
According to the references [33][34][35][36], considering the limit value of the migration of fine particles [7,[37][38][39]. The migration equation of fine particles is as follows: where: φ max is the maximum porosity, which refers to the limit value of porosity change in porous media caused by the migration of fine particles; C cr is the critical volume fraction of fine particles when the migration reaches equilibrium; |v l | is the modulus of the velocity vector of the fluid (m/s); ζ is the erosion coefficient of fine particles (m −1 ).

Auxiliary Equation
(1) Porosity evolution equation The change in pore volume is mainly caused by the erosion and migration of fine particles and the deformation of solid skeleton. For the erosion and migration of fine particles, the rate of change in pore volume is equal to the volume of eroded fine particles per unit time, then: Appl. Sci. 2020, 10, 7539 Which can be simplified as: For pore volume change caused by deformation of solid skeleton, the porosity equation proposed in reference [40] is as follows: where, φ 0 is the initial porosity. Equation (23) can be written as: Then, based on the erosion and migration of fine particles and the deformation of solid skeleton, the porosity evolution equation can be expressed as: (2) Relationship between permeability and porosity The Kozeny-Carman equation was used to describe the relationship between permeability and porosity [41][42][43]: where: k 0 is the initial permeability (m 2 ). Equations (3)~(5), (11), (19), (20), (25) and (26), constitute the two-phase fluid-solid coupling mathematical model considering the effects of fine particles migration.

Boundary Conditions and Initial Conditions
The coupling mathematical model mentioned above is highly nonlinear, boundary conditions and initial conditions need to be given in the process of solution. The pore water pressure, the volume fraction of fine particles, and the displacement of solid skeleton are taken as basic variables.
Boundary conditions: the boundary conditions usually have the first type, also known as Dirichlet boundary; and the second type, also known as Neumann boundary.
The first boundary conditions can be expressed as: The second boundary conditions can be expressed as: q l = −v ls · n on Γ q l ; q C = −Cv ps · n on Γ C ; σ ij · n = t on Γ σ ij . Initial conditions: p l = p 0 l , C = C 0 , u s = u 0 s at t = 0.

The Solution of the Mathematical Model
The above mathematical model was solved by using the coupling solid mechanics module and porous medium multiphase flow module in the software COMSOL Multiphysics ® . Taking the volume fraction of fine particles (C), the pore pressure of the fluid (p l ), and the displacement of the solid skeleton (u s ) as the basic unknowns, the other unknowns such as the velocity vector of the fluid v l , the velocity vector of the fine particles (v p ), the mass rate of fine particles migration ( . m), the porosity (φ), and the permeability (k 0 )can be calculated by the basic unknowns.

Validation of the Coupling Model
In this section, the software COMSOL Multiphysics ® is used to solve the mathematical model. Moreover, the solution results were compared with those in reference [33] to verify the accuracy of the coupling model in this paper.

The Establishment of Numerical Model
Reference [33] studied the problem of sand production in the wellbore during oil extraction. A two-dimensional calculation model is adopted. The radius of the wellbore (inner boundary) is 0.1 m, and the outer boundary is 5 m. The calculation model is shown in Figure 3.  The standard atmospheric pressure value is taken as the calculation reference value. The initial porosity of the model was 0.25, and the initial volume fraction of fine particles was 0.001. The boundary conditions adopted by the model are as follows: the fluid pressure at inner boundary is 5 MPa, and the fluid pressure at outer boundary is 8 MPa. The inner boundary is the outflow boundary for fine particles, and the external stress acting on the outer boundary is 20 MPa. The physical parameters are shown in Table 1.  The standard atmospheric pressure value is taken as the calculation reference value. The initial porosity of the model was 0.25, and the initial volume fraction of fine particles was 0.001. The boundary conditions adopted by the model are as follows: the fluid pressure at inner boundary is 5 MPa, and the fluid pressure at outer boundary is 8 MPa. The inner boundary is the outflow boundary for fine particles, and the external stress acting on the outer boundary is 20 MPa. The physical parameters are shown in Table 1.

Verification and Analysis of Calculation Results
The evolution curves of volume fraction of fine particles, porosity, velocity of fluid, and radial displacement at the inner boundary calculated according to the mathematical model in this paper are shown in Figures 4-7 respectively, and the calculation results in this paper are compared with those in reference [33]. A comprehensive analysis of Figures 4-7 shows that the evolution curves of volume fraction of fine particles, porosity, velocity of fluid, and radial displacement calculated by the mathematical model established in this paper are consistent to the results in reference [33]. Because of the difference between the porosity evolution equation, the permeability equation in this paper, and those in reference [33], there is a certain error between the calculation results and those in reference [33]. However, the error is within the acceptable range. This also verifies the correctness of the model in this paper.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 24 of the difference between the porosity evolution equation, the permeability equation in this paper, and those in reference [33], there is a certain error between the calculation results and those in reference [33]. However, the error is within the acceptable range. This also verifies the correctness of the model in this paper.   id v l (m/s) Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 24 of the difference between the porosity evolution equation, the permeability equation in this paper, and those in reference [33], there is a certain error between the calculation results and those in reference [33]. However, the error is within the acceptable range. This also verifies the correctness of the model in this paper.

The Establishment of Numerical Model
According to "Code for Design of Heavy Haul Railway (TB10625-2017)" [44] in China and referring to a typical section structure of the Shuo-huang Railway subgrade, a two-dimensional calculation model is established as shown in Figure 8. The structure of the subgrade from top to bottom is as follows: surface layer of the subgrade (the thickness is 0.6 m), bottom layer of the subgrade (the thickness is 1.9 m). The slope gradient of the subgrade is 1:1.5. The model sets four observation points, as shown in Figure 8. According to references [45][46][47][48], we use the free mesh triangular elements and an ultra-fine mesh (the number of mesh is 922, the minimum mass for the mesh is 0.55 and the average mass for the mesh is 0.89) is set in this research. The mesh of the calculation model is shown in Figure 9a, and the influence of mesh number on the velocity of fluid and the influence of mesh number on the vertical displacement is shown in Figure 9b. The results indicate that, the calculation results of models with the mesh number of 138, 303, and 922 are different. When the mesh number is increased to 922, the calculation results are little different from those of

The Establishment of Numerical Model
According to "Code for Design of Heavy Haul Railway (TB10625-2017)" [44] in China and referring to a typical section structure of the Shuo-huang Railway subgrade, a two-dimensional calculation model is established as shown in Figure 8. The structure of the subgrade from top to bottom is as follows: surface layer of the subgrade (the thickness is 0.6 m), bottom layer of the subgrade (the thickness is 1.9 m). The slope gradient of the subgrade is 1:1.5. The model sets four observation points, as shown in Figure 8. According to references [45][46][47][48], we use the free mesh triangular elements and an ultra-fine mesh (the number of mesh is 922, the minimum mass for the mesh is 0.55 and the average mass for the mesh is 0.89) is set in this research. The mesh of the calculation model is shown in Figure 9a, and the influence of mesh number on the velocity of fluid and the influence of mesh number on the vertical displacement is shown in Figure 9b. The results indicate that, the calculation results of models with the mesh number of 138, 303, and 922 are different. When the mesh number is increased to 922, the calculation results are little different from those of model with the mesh number 3566. Therefore, the ultra-fine mesh in this paper is accepted.
subgrade (the thickness is 1.9 m). The slope gradient of the subgrade is 1:1.5. The model sets four observation points, as shown in Figure 8. According to references [45][46][47][48], we use the free mesh triangular elements and an ultra-fine mesh (the number of mesh is 922, the minimum mass for the mesh is 0.55 and the average mass for the mesh is 0.89) is set in this research. The mesh of the calculation model is shown in Figure 9a, and the influence of mesh number on the velocity of fluid and the influence of mesh number on the vertical displacement is shown in Figure 9b. The results indicate that, the calculation results of models with the mesh number of 138, 303, and 922 are different. When the mesh number is increased to 922, the calculation results are little different from those of model with the mesh number 3566. Therefore, the ultra-fine mesh in this paper is accepted.

Mechanical Parameters of Subgrade Filler
The filler of the subgrade is made of graded crushed stone. In order to analyze the influence of volume fraction of fine particles on the fluid-solid two-phase coupling characteristics of the filler, three different volume fractions of fine particles are selected, and the gradation of the graded crushed stone filler meet the requirements of the specification [49]. The physical parameters are shown in Table 2.

Mechanical Parameters of Subgrade Filler
The filler of the subgrade is made of graded crushed stone. In order to analyze the influence of volume fraction of fine particles on the fluid-solid two-phase coupling characteristics of the filler, three different volume fractions of fine particles are selected, and the gradation of the graded crushed stone filler meet the requirements of the specification [49]. The physical parameters are shown in Table 2.

Boundary Conditions and Initial Conditions
The initial conditions for the calculation model conclude that the pore water pressure is the hydrostatic pressure; the initial volume fractions of fine particles are 0.01, 0.03, and 0.07, respectively; the initial porosity of the subgrade is 0.3; and the initial displacement is 0.
It is assumed that the subgrade in the calculated area is saturated, and the dissipation path of pore water pressure is from bottom to top. The seepage boundary conditions for the calculation model are as follows: the lower boundary AD is the boundary of fluid pressure, and the fluid pressure at the boundary AD is set as 40 kPa according to the reference [1][2][3]. The fluid pressure at upper boundary BC is 0 kPa. The slopes are usually protected, so AB and CD are assumed to be impermeable boundaries. The boundary conditions corresponding to the migration of fine particles are as follows: the boundary AD is the inflow boundary for fine particles carried by the fluid, and the volume fractions of fine particles at the boundary AD are 0.01, 0.03, and 0.07, respectively. The upper boundary BC is the outflow boundary of fine particles. AB and CD are no flux boundary. The displacement boundary conditions are as follows: the boundary AD is set as a fixed boundary, the remaining boundaries are free boundaries.  Figure 11 shows the evolution curves of porosity in saturated zone of railway subgrade, in which subfigure (b) is the detailed curves for the curves in the box of subfigure (a). It can be seen from Figure  11 that:

The Variation of Porosity
(1) The variation trend of porosity at the four observation points is same. First, it increases slowly, then it increases rapidly to a certain value (about 0.5, i.e., the maximum porosity), and then it is basically stable. (2) For observation points 1 and 2, the time for the porosity to reach the maximum porosity is same; for observation point 3, the time for the porosity to reach the maximum porosity is slightly later than that of observation points 1 and 2; while for observation point 4, it has the latest time to reach the maximum porosity.  Figure 11 shows the evolution curves of porosity in saturated zone of railway subgrade, in which subfigure (b) is the detailed curves for the curves in the box of subfigure (a). It can be seen from Figure 11 that:

The Variation of Porosity
(1) The variation trend of porosity at the four observation points is same. First, it increases slowly, then it increases rapidly to a certain value (about 0.5, i.e., the maximum porosity), and then it is basically stable. (2) For observation points 1 and 2, the time for the porosity to reach the maximum porosity is same; for observation point 3, the time for the porosity to reach the maximum porosity is slightly later than that of observation points 1 and 2; while for observation point 4, it has the latest time to reach the maximum porosity.  Figure 11 shows the evolution curves of porosity in saturated zone of railway subgrade, in which subfigure (b) is the detailed curves for the curves in the box of subfigure (a). It can be seen from Figure  11 that:

The Variation of Porosity
(1) The variation trend of porosity at the four observation points is same. First, it increases slowly, then it increases rapidly to a certain value (about 0.5, i.e., the maximum porosity), and then it is basically stable. (2) For observation points 1 and 2, the time for the porosity to reach the maximum porosity is same; for observation point 3, the time for the porosity to reach the maximum porosity is slightly later than that of observation points 1 and 2; while for observation point 4, it has the latest time to reach the maximum porosity.  Figure 12 shows the distribution cloud map of pore water pressure in saturated zone of railway subgrade, and Figure 13 shows the evolution curve of pore water pressure at the observation points 3. From these two pictures we can see that:

The Variation of Pore Water Pressure
(1) The pore water pressure at the lower boundary of the subgrade is 65 kPa, and the pore water pressure at the upper boundary of the subgrade is 0 kPa. The pore water pressure is evenly distributed from the bottom of the subgrade to the top of the subgrade. (2) With the increase of time, the pore water pressure at point 3 shows a trend of first decreasing, then increasing, and then decreasing, and finally remains stable, but the range of its increase or decrease is small. This reason is that with the erosion and migration of fine particles in the subgrade, the porosity and permeability of each point inside the subgrade are different (see Figures 10 and 14 for details), thus affecting the distribution of pore water pressure. Therefore, the pore water pressure at the same point is not constant.   Figure 12 shows the distribution cloud map of pore water pressure in saturated zone of railway subgrade, and Figure 13 shows the evolution curve of pore water pressure at the observation points 3. From these two pictures we can see that:

The Variation of Pore Water Pressure
(1) The pore water pressure at the lower boundary of the subgrade is 65 kPa, and the pore water pressure at the upper boundary of the subgrade is 0 kPa. The pore water pressure is evenly distributed from the bottom of the subgrade to the top of the subgrade. (2) With the increase of time, the pore water pressure at point 3 shows a trend of first decreasing, then increasing, and then decreasing, and finally remains stable, but the range of its increase or decrease is small. This reason is that with the erosion and migration of fine particles in the subgrade, the porosity and permeability of each point inside the subgrade are different (see Figures 10 and 14 for details), thus affecting the distribution of pore water pressure. Therefore, the pore water pressure at the same point is not constant.   Figure 12 shows the distribution cloud map of pore water pressure in saturated zone of railway subgrade, and Figure 13 shows the evolution curve of pore water pressure at the observation points 3. From these two pictures we can see that:

The Variation of Pore Water Pressure
(1) The pore water pressure at the lower boundary of the subgrade is 65 kPa, and the pore water pressure at the upper boundary of the subgrade is 0 kPa. The pore water pressure is evenly distributed from the bottom of the subgrade to the top of the subgrade. (2) With the increase of time, the pore water pressure at point 3 shows a trend of first decreasing, then increasing, and then decreasing, and finally remains stable, but the range of its increase or decrease is small. This reason is that with the erosion and migration of fine particles in the subgrade, the porosity and permeability of each point inside the subgrade are different (see Figures 10 and 14 for details), thus affecting the distribution of pore water pressure. Therefore, the pore water pressure at the same point is not constant.  (1) The variation trend of fluid velocity is basically the same as that of fine particles migration velocity, showing a trend of basically unchanged at first, then rapidly increasing, and finally tending to be stable. (2) The velocity of fluid is always greater than the velocity of fine particles, and the velocity difference between these two is becoming more and more obvious with the increase of fluid velocity. The reason is that with the increase of the fluid velocity, the drag coefficient decreases, and the drag force of the fine particles is reduced, which hinders the movement of the fine particles to some extent. Figure 15 shows the relation curves between fluid velocity and porosity at different observation points. It can be seen from Figure 15 that the difference between the fluid velocity and porosity at different observation points is small. The fluid velocity increases with the increase of porosity, and shows an exponential change.  showing a trend of basically unchanged at first, then rapidly increasing, and finally tending to be stable. (2) The velocity of fluid is always greater than the velocity of fine particles, and the velocity difference between these two is becoming more and more obvious with the increase of fluid velocity. The reason is that with the increase of the fluid velocity, the drag coefficient decreases, and the drag force of the fine particles is reduced, which hinders the movement of the fine particles to some extent. Figure 15 shows the relation curves between fluid velocity and porosity at different observation points. It can be seen from Figure 15 that the difference between the fluid velocity and porosity at different observation points is small. The fluid velocity increases with the increase of porosity, and shows an exponential change.  (1) The variation trend of fluid velocity is basically the same as that of fine particles migration velocity, showing a trend of basically unchanged at first, then rapidly increasing, and finally tending to be stable. (2) The velocity of fluid is always greater than the velocity of fine particles, and the velocity difference between these two is becoming more and more obvious with the increase of fluid velocity. The reason is that with the increase of the fluid velocity, the drag coefficient decreases, and the drag force of the fine particles is reduced, which hinders the movement of the fine particles to some extent. Figure 15 shows the relation curves between fluid velocity and porosity at different observation points. It can be seen from Figure 15 that the difference between the fluid velocity and porosity at different observation points is small. The fluid velocity increases with the increase of porosity, and shows an exponential change.   Figure 16 shows the evolution curves of permeability in saturated zone of railway subgrade. It can be seen from Figure 16 that the permeability increases rapidly to the maximum after a relatively slow growth, and finally remains stable. For observation points 1~3, the maximum permeability is about 4.9 × 10 −11 m 2 , which is an increase of about 6.3 times higher than the initial permeability.  Figure 17 shows the evolution curves of the volume fractions of fine particles with different initial fine particles volume fraction at point 2. From Figure 17, we know that the larger the initial fine particles volume fraction is, the larger the maximum value of fine particles volume fraction is, and the shorter the time needed to reach the maximum value. Combined with Equation (20), the larger the initial fine particles volume fraction is, the larger the mass rate of fine particles migration will be; therefore, the shorter the time needed to reach the maximum fine particles volume fractions will be. It can be concluded that the fine particles content in the filler of subgrade should be strictly controlled, and the fine particles content should be reduced as far as possible on the premise that the filler grade meets the requirements of the specification in practical engineering.  Figure 16 shows the evolution curves of permeability in saturated zone of railway subgrade. It can be seen from Figure 16 that the permeability increases rapidly to the maximum after a relatively slow growth, and finally remains stable. For observation points 1~3, the maximum permeability is about 4.9 × 10 −11 m 2 , which is an increase of about 6.3 times higher than the initial permeability.  Figure 16 shows the evolution curves of permeability in saturated zone of railway subgrade. It can be seen from Figure 16 that the permeability increases rapidly to the maximum after a relatively slow growth, and finally remains stable. For observation points 1~3, the maximum permeability is about 4.9 × 10 −11 m 2 , which is an increase of about 6.3 times higher than the initial permeability.  Figure 17 shows the evolution curves of the volume fractions of fine particles with different initial fine particles volume fraction at point 2. From Figure 17, we know that the larger the initial fine particles volume fraction is, the larger the maximum value of fine particles volume fraction is, and the shorter the time needed to reach the maximum value. Combined with Equation (20), the larger the initial fine particles volume fraction is, the larger the mass rate of fine particles migration will be; therefore, the shorter the time needed to reach the maximum fine particles volume fractions will be. It can be concluded that the fine particles content in the filler of subgrade should be strictly controlled, and the fine particles content should be reduced as far as possible on the premise that the filler grade meets the requirements of the specification in practical engineering.  Figure 17 shows the evolution curves of the volume fractions of fine particles with different initial fine particles volume fraction at point 2. From Figure 17, we know that the larger the initial fine particles volume fraction is, the larger the maximum value of fine particles volume fraction is, and the shorter the time needed to reach the maximum value. Combined with Equation (20), the larger the initial fine particles volume fraction is, the larger the mass rate of fine particles migration will be; therefore, the shorter the time needed to reach the maximum fine particles volume fractions will be. It can be concluded that the fine particles content in the filler of subgrade should be strictly controlled, and the fine particles content should be reduced as far as possible on the premise that the filler grade meets the requirements of the specification in practical engineering. Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 24  Figure 18 shows the contours of vertical effective stress at the steady state of fluid seepage. Figure 19 shows the evolution curves of vertical effective stress at different observation points in saturated zone of railway subgrade. It can be seen from Figure 18 that the vertical effective stress is distributed symmetrically on the central axis of the subgrade. With the increase of depth, the vertical effective stress decreases. It can be seen from Figure 19 that vertical effective stress at points 2 and points 3 increases first and then decreases slightly with the increase of time, and finally tends to be stable. In combination with Figure 13, the pore water pressure increases at a time of about 1000-2000 min, so the effective stress decreases.   Figure 18 shows the contours of vertical effective stress at the steady state of fluid seepage. Figure 19 shows the evolution curves of vertical effective stress at different observation points in saturated zone of railway subgrade. It can be seen from Figure 18 that the vertical effective stress is distributed symmetrically on the central axis of the subgrade. With the increase of depth, the vertical effective stress decreases. It can be seen from Figure 19 that vertical effective stress at points 2 and points 3 increases first and then decreases slightly with the increase of time, and finally tends to be stable. In combination with Figure 13, the pore water pressure increases at a time of about 1000-2000 min, so the effective stress decreases.  Figure 18 shows the contours of vertical effective stress at the steady state of fluid seepage. Figure 19 shows the evolution curves of vertical effective stress at different observation points in saturated zone of railway subgrade. It can be seen from Figure 18 that the vertical effective stress is distributed symmetrically on the central axis of the subgrade. With the increase of depth, the vertical effective stress decreases. It can be seen from Figure 19 that vertical effective stress at points 2 and points 3 increases first and then decreases slightly with the increase of time, and finally tends to be stable. In combination with Figure 13, the pore water pressure increases at a time of about 1000-2000 min, so the effective stress decreases.   Figure 20 shows the contours of vertical displacement at the steady state of fluid seepage. We can see that the vertical displacement is distributed symmetrically on the central axis of the subgrade. At the same depth, the vertical displacement at the central axis of the subgrade is slightly larger than the vertical displacement at the edge of the subgrade. Figure 21 shows the evolution curves of vertical displacement at different observation points in saturated zone of railway subgrade. It can be seen from Figure 21 Figure 20 shows the contours of vertical displacement at the steady state of fluid seepage. We can see that the vertical displacement is distributed symmetrically on the central axis of the subgrade. At the same depth, the vertical displacement at the central axis of the subgrade is slightly larger than the vertical displacement at the edge of the subgrade. Figure 21 shows the evolution curves of vertical displacement at different observation points in saturated zone of railway subgrade. It can be seen from Figure 21 Figure 20 shows the contours of vertical displacement at the steady state of fluid seepage. We can see that the vertical displacement is distributed symmetrically on the central axis of the subgrade. At the same depth, the vertical displacement at the central axis of the subgrade is slightly larger than the vertical displacement at the edge of the subgrade. Figure 21 shows the evolution curves of vertical displacement at different observation points in saturated zone of railway subgrade. It can be seen from Figure 21 Figure 22a shows the relation curves of porosity change and accumulative volume fraction of fine particles and (b) shows the relation curves of porosity change and vertical displacement. The positive value of porosity change represents the increase of porosity, and the negative value of porosity change represents the decrease of porosity. It can be seen from Figure 22 that the accumulative volume fraction of fine particles has a significant impact on the change of porosity: with the increase of the accumulation volume fraction of fine particles, porosity change increases gradually, when the accumulation volume fraction of fine particles reaches about 10, the porosity change reaches a plateau value. The vertical displacement has little influence on the change of porosity. With the change of vertical displacement, the porosity first decreases slightly and then increases slightly, and finally remains stable.  Figure 22a shows the relation curves of porosity change and accumulative volume fraction of fine particles and (b) shows the relation curves of porosity change and vertical displacement. The positive value of porosity change represents the increase of porosity, and the negative value of porosity change represents the decrease of porosity. It can be seen from Figure 22 that the accumulative volume fraction of fine particles has a significant impact on the change of porosity: with the increase of the accumulation volume fraction of fine particles, porosity change increases gradually, when the accumulation volume fraction of fine particles reaches about 10, the porosity change reaches a plateau value. The vertical displacement has little influence on the change of porosity. With the change of vertical displacement, the porosity first decreases slightly and then increases slightly, and finally remains stable.

Conclusions
Based on the theory of two-phase seepage, the theory of porous media seepage, and the principle of effective stress in porous media, a two-phase fluid-solid coupling mathematical model in saturated zone of subgrade considering the effects of fine particles migration is established. The mathematical model is numerically calculated with the software COMSOL Multiphysics ® . The two-phase seepage characteristics in saturated zone of the subgrade and the deformation characteristics of the solid skeleton are studied, and the following conclusions are obtained: (1) Under the action of fluid seepage, the erosion and migration of the fine particles in saturated zone of subgrade gradually occurred, and the volume fraction of fine particles increased with the increase of time, and as the time goes on, only skeleton particles and particles with large sizes are left in the subgrade, the volume fraction of fine particles gradually decreased to the initial value.

Conclusions
Based on the theory of two-phase seepage, the theory of porous media seepage, and the principle of effective stress in porous media, a two-phase fluid-solid coupling mathematical model in saturated zone of subgrade considering the effects of fine particles migration is established. The mathematical model is numerically calculated with the software COMSOL Multiphysics ® . The two-phase seepage characteristics in saturated zone of the subgrade and the deformation characteristics of the solid skeleton are studied, and the following conclusions are obtained: (1) Under the action of fluid seepage, the erosion and migration of the fine particles in saturated zone of subgrade gradually occurred, and the volume fraction of fine particles increased with the increase of time, and as the time goes on, only skeleton particles and particles with large sizes are left in the subgrade, the volume fraction of fine particles gradually decreased to the initial value. The constant erosion of fine particles inside the subgrade and the constant migration of fine particles with fluid from the lower part of the subgrade to the upper part of the subgrade caused that the volume fraction of fine particles in the upper part of the subgrade is larger than that in the lower part of the subgrade at the same time.
(2) With the increase of time, the porosity, the velocity of fluid, the velocity of fine particles, and the permeability of the subgrade show the trend of first increasing and then stabilizing; the pore water pressure in the subgrade has no significant changes with time. (3) The vertical displacement is distributed symmetrically on the central axis of the subgrade. With the increase of time, the vertical displacement increases first and then decreases slightly and finally tends to be stable. (4) For the filler with a larger initial volume fraction of fine particles, the maximum value of the volume fraction of fine particles in the subgrade caused by fluid seepage is larger, and the time required to reach the maximum value is shorter. It can be concluded that in actual engineering, the volume fraction of fine particles in the subgrade filler should be minimized on the premise so that the filler gradation meets the requirements of the specification.
The main purpose of this research is to explore the characteristics of the two-phase seepage and the deformation of solid skeleton in saturated zone of the subgrade considering the effects of fine particles migration, in order to explain the internal mechanisms of the mud pumping disease. Therefore, the dynamic characteristics of the subgrade filler and the fluid-solid two-phase seepage characteristics considering the effects of fine particles migration and sedimentation under train loads have more practical engineering significance, which will also be the focus of the next research.