\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[greek,ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{Study on Particle-Size Process on Internal Erosion of Grap-Graded
Soil--Rock Mixtures of Different Fine Contents}



\author[1]{Zhilin Cao}%
\author[1]{Weichen Sun}%
\author[1]{qiang xie}%
\author[1]{Zhihui Wu}%
\author[2]{Xiang Fu}%
\author[1]{Dalang Tian}%
\affil[1]{Chongqing University}%
\affil[2]{Chongqing Jiaotong University}%


\vspace{-1em}



  
  \date{March 12, 2025}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Soil--rock mixtures are widely encountered in geotechnical engineering
projects. The instability and failure mechanism of grap-graded
soil--rock mixtures under rainfall conditions has always been the focus
of geological disaster research. To deeply explore the mechanism of
seepage deformation of soil--rock mixtures, an indoor physical
permeability test that considers soil--rock mixtures with different fine
contents was conducted, and a particle-scale numerical simulation test
of the permeability evolution was carried out using the coupling model
of PFC3D and ABAQUS. The test results showed that the spatial
distribution of fine particle loss along the height direction could be
divided into three areas: top loss, middle uniform, and bottom loss
area. The ``island'' effect of coarse particles, which is caused by
excessive fine content and makes the fine particles bear more load, was
eliminated with the loss of fine particles. In this preset working
condition of coarse and fine particle diameters, setting FC to 35\% may
be the best way to fill the voids between the coarse particles. Particle
migration leads to a change in the load-bearing skeleton structure,
thereby causing seepage deformation. Therefore, the particle-scale
numerical test method can better reproduce the seepage deformation
process of grap-graded soil--rock mixtures.%
\end{abstract}%



\sloppy


\textbf{Study on Particle-Size Process on Internal Erosion of
Grap-Graded Soil--Rock Mixtures of Different Fine Contents}

\textbf{Zhilin Cao\textsuperscript{1}, Weichen Sun\textsuperscript{1},
Qiang Xie\textsuperscript{1*}, Zhihui Wu\textsuperscript{1}, Xiang
Fu\textsuperscript{2}, Dalang Tian\textsuperscript{1}}

\textsuperscript{1}School of Civil Engineering, Chongqing University,
Chongqing, China\textsuperscript{2}College of River and Ocean
Engineering, Chongqing Jiaotong University, Chongqing, China

\textbf{Corresponding Author} :

address:School of Civil Engineering, Chongqing University, Chongqing,
China

Email: xieqiang2000@163.com

\textbf{Abstract}

Soil--rock mixtures are widely encountered in geotechnical engineering
projects. The instability and failure mechanism of grap-graded
soil--rock mixtures under rainfall conditions has always been the focus
of geological disaster research. To deeply explore the mechanism of
seepage deformation of soil--rock mixtures, an indoor physical
permeability test that considers soil--rock mixtures with different fine
contents was conducted, and a particle-scale numerical simulation test
of the permeability evolution was carried out using the coupling model
of PFC3D and ABAQUS. The test results showed that the spatial
distribution of fine particle loss along the height direction could be
divided into three areas: top loss, middle uniform, and bottom loss
area. The ``island'' effect of coarse particles, which is caused by
excessive fine content and makes the fine particles bear more load, was
eliminated with the loss of fine particles. In this preset working
condition of coarse and fine particle diameters, setting FC to 35\% may
be the best way to fill the voids between the coarse particles. Particle
migration leads to a change in the load-bearing skeleton structure,
thereby causing seepage deformation. Therefore, the particle-scale
numerical test method can better reproduce the seepage deformation
process of grap-graded soil--rock mixtures.

\textbf{Keywords}

Soil--rock mixtures, internal erosion, fines content, solid-fluid
interaction, particle skeleton structure

\subsection*{INTRODUCTION}

{\label{introduction}}

Soil--rock mixtures consist of rock blocks within soils composed of
sand, silt, and clay. The geological origins of soil--rock mixtures vary
between landslide, colluvial, alluvial--diluvial, and glacial deposits
(Xu et al., 2007), which are widely distributed in the Three Gorges
Reservoir area, China (Zhang et al., 2015; Shi et al., 2016). Under the
action of rainfall infiltration, landslides, collapses, and other
geological disasters can easily occur due to the migration of fine
particles in gap-graded soil--rock mixtures (Grosta et al., 1999; Zhang
et al., 2011; Johan et al., 2013; Lei et al., 2017; Fan et al., 2020).
Therefore, the study of seepage-induced fine particle migration through
gap-graded soils is crucial (Cheng et al., 2018; Tian et al., 2020).

Soil--rock mixtures are granular materials that can be constituted in
many different ways with different types of intergrain contacts
(Thevanayagam et al., 2002). The microstructure and mechanical behaviour
of granular materials largely depend on the particles of the constituent
materials and the particle size distributions (PSDs). Particles of
different sizes exhibit a high degree of force heterogeneity. Skempton
et al. (1994) found that fine content (FC) in gap-graded soils affects
soil strength. Ouyang et al. (2015) presented experimental
investigations on the influence of initial FC on the fabric of soils
subjected to internal erosion. Lopez et al. (2016) proposed a case study
of a soil structure based on previous empirical observations, which
conceptually explained the effects of different components (coarse and
fine particles) on the mixture behaviour. Tian et al. (2020) studied the
critical FC of natural accumulation soil, beyond which the internal
stability of the mix was distinctly altered. Based on the fact that fine
particles in soils gradually migrate downwards and lead to subsequent
slope failure during rainfall events, Yin et al. (2021) developed a 3D
discrete solid--fluid sequentially coupled model to analyse changes in
pore structure caused by fine particle migration. However, there are few
studies on the micro process and evolution characteristics of fine
particle migration of different FCs and the visual description of the
infiltration process phenomena.

Owing to advancements in computer technology, numerical simulation
methods have been employed to address this limitation and investigate
the particle migration process with high precision (Guo and Yu, 2017;
Peters et al., 2019; Wang et al., 2020). For the study of gap-graded
soils, particle-size numerical test using the CFD--DEM method has become
a research trend (Zou et al., 2020; Xiong et al., 2021), as it can
enable a deeper understanding of the fine particle migration process and
the coupling mechanism between particles and seepage flow. Zou et al.
(2013) studied the particle transport mechanism in a base soil--filter
system using the CFD--DEM method. Yang et al. (2019) improved the
semi-resolved CFD--DEM model for seepage-induced fine particle
migration. Cheng et al. (2021) revealed that unresolved CFD--DEM models
predict unreasonably high critical hydraulic gradients for suffusion in
gap-graded soil columns packed under gravity.

The occurrence of internal erosion of grap-graded soil--rock mixtures is
the result of the interaction of water pressure, particle migration, and
soil stress state. In this study, a fluid--solid coupled model of
internal erosion through the coupling of CFD and DEM was developed to
study the coupling effect of water pressure, particle migration, and
soil stress state. Moreover, a permeability test for the stability of
the gap-graded soil--rock mixtures considering particle loss was
conducted to analyse the effect of different FCs, and the numerical
simulation results were verified. This study focused on the
particle-scale process of the seepage deformation of grap-graded
soil--rock mixtures and dynamic changes in load-bearing skeletons with
different FCs. This study is helpful in exploring the deformation and
failure law of soil--rock mixtures under the condition of rainfall
infiltration.

\subsection*{COUPLED CFD--DEM MODEL}

{\label{coupled-cfddem-model}}

The coupled CFD--DEM method used in this study included three types of
formulation: DEM, CFD, and CFD--DEM coupling formulations. PFC3D was
used to solve the particle system at the particle scale, and Newton's
law of motion was used to control the particle motion. Moreover, ABAQUS
with the porous media seepage theory, Darcy flow model, was used as the
CFD solver. The particle--fluid interaction, including the drag and
pressure gradient forces, was calculated using the ABAQUS--PFC3D program
developed by us. The entire fluid--solid interaction process is
illustrated in Fig. 1.

\subsubsection*{2.1 Governing Equations of the Particle
Phase}

{\label{governing-equations-of-the-particle-phase}}

In the fluid--solid coupling method, the behaviour of the particle phase
is controlled by Newton's law of motion. In the PFC3D program, the
governing equations of the particle phase are given by Eqs. (1) and (2)
(Itasca, 2018),\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

(1)\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

(2)

where u is the particle velocity, m is the particle mass,
f\textsubscript{fluid} is the total force applied by the fluid on the
particle, f\textsubscript{mesh} is the sum of additional forces
(externally applied forces and contact forces) acting on the particle, g
is the acceleration due to gravity, \selectlanguage{greek}ω \selectlanguage{english}is the particle angular velocity,
I is the moment of inertia, and M is the moment acting on the particle.

\subsubsection*{2.2 Governing Equations of the Fluid
Phase}

{\label{governing-equations-of-the-fluid-phase}}

In the CFD--DEM method, the Darcy fluid model has been proven to be an
effective and efficient method for predicting fluid behaviour using an
appropriate permeability (Goodarzi et al., 2014; Cui et al., 2016; Tang
et al., 2017). The Darcy fluid model was adopted for the fluid flow
simulation in this study. The governing equations of the transient fluid
flow through a saturated porous media based on the Darcy model and the
fluid phase are assumed to be incompressible.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textsubscript{(3)}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textsubscript{(4)}

where v is the fluid velocity, k is the matrix permeability,
\selectlanguage{greek}μ\selectlanguage{english}\textsubscript{f} is the fluid viscosity, p is the fluid pressure,
\selectlanguage{greek}ρ\selectlanguage{english}\textsubscript{f} is the fluid density, g is the magnitude of
gravitational acceleration, and Z is the position height. Given the
boundary conditions, the fluid pressure field can be calculated using
Eq. (4), and the corresponding water velocity can be determined using
Eq. (3).

\subsubsection*{2.3 Fluid--Particle Interaction
Force}

{\label{fluidparticle-interaction-force}}

The fluid--solid interaction follows the numerical analysis system of
the coarse-grid approach (Tsuji et al., 1993), as shown in Fig. 2. In
the coarse-grid method, equations describing the fluid flow are solved
numerically on a set of elements that are larger than the particles. The
force applied by the fluid to the particles (the fluid--particle
interaction force), f\textsubscript{fluid} , comprises two parts: a drag
force and a force due to the fluid pressure gradient, which are given in
Eq. (5).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

(5)

The drag force exerted by the fluid on the particles is defined
individually for each particle. Depending on how the porosity is
calculated, a particle may intersect with more than one fluid element.
The fluid--solid interaction force is always applied at the particle
centroid, and no rotational moment is applied to the particle. The drag
force, f\textsubscript{drag} , is defined as:\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

(6)

where f\textsubscript{0} is a single-particle drag force given by Eq.
(7) and \selectlanguage{greek}ε \selectlanguage{english}is the porosity of the fluid element where the particle
resides. The \selectlanguage{greek}ε\selectlanguage{english}\selectlanguage{greek}\textsuperscript{-χ} \selectlanguage{english}term is an empirical factor that
accounts for local porosity. This correction term makes the force
applicable to both high-porosity and low-porosity systems and for a
large range of Reynolds numbers (Difelice et al., 1994; Xu et al.,
1997).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

(7)

where C\textsubscript{d} is the drag
coefficient,\selectlanguage{greek}\emph{ρ\textsubscript{f}}\selectlanguage{english} is the fluid density, r is the
particle radius, v is the fluid velocity, and u is the particle
velocity.

The drag coefficient is defined as:\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}

\end{center}
\end{figure}

(8)

where R\textsubscript{ep} is the particle Reynolds number.

The empirical coefficient \selectlanguage{greek}χ \selectlanguage{english}is defined as:\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image9/image9}

\end{center}
\end{figure}

(9)

The particle Reynolds number is given as:\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image10/image10}

\end{center}
\end{figure}

(10)

where \selectlanguage{greek}\emph{μ\textsubscript{f}}\selectlanguage{english} is the fluid viscosity.

The Kozeny--Carman relation (Kozeny et al., 1927; Carman 1939) was used
to estimate the permeability. Two-way coupling was accomplished by
updating the permeability coefficient information for each fluid element
and updating the fluid velocity field of each discrete element. The
permeability of each fluid element was calculated from the porosity of
the PFC3D model, and the calculation is given by Eq. (11),\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image11/image11}

\end{center}
\end{figure}

(11)

where k is the permeability, r is the radius of the PFC particles, c is
the K--C constant, which was 5 in this study, and \selectlanguage{greek}ε\selectlanguage{english}\textsuperscript{i}
is the porosity of each fluid element. The permeability coefficient was
taken as the upper limit of porosity, which was 70\%. When the porosity
of the element exceeded 70\%, most particles in the fluid element were
washed away, and the permeability coefficient was taken as a fixed
value.

\subsection*{PHYSICAL TEST}

{\label{physical-test}}

A vertical erosion test device was developed for this study. The device
can simulate the action of a constant water head and load on a sample in
the saturated state to ensure the formation of a dynamic seepage network
in soil--rock mixtures.

\subsubsection*{3.1 Device Design}

{\label{device-design}}

The set of devices included a sample saturation device, an erosion
simulation device with a constant head, a particle collection device,
and data collection devices. The erosion simulation device was divided
into five layers, which was convenient for analysing the changes in FC
at different heights after the test, as shown in Fig. 3.

The sample saturation device was a transparent plexiglass water tank.
The water tank was composed of five transparent plexiglass plates with a
thickness of 1 cm. The erosion simulation device was composed of five
transparent plexiglass cylinders, which could be easily disassembled and
assembled. The fifth plexiglass cylinderwas a constant head control part
with a height of 25 cm, and an overflow drain was installed 5 cm away
from the top of the fifth cylinder to keep the head height constant at
10 cm (the global hydraulic gradient, i , was 0.2). The five transparent
plexiglass cylinders were installed by a flange connection, and they
were sealed and impervious during the test.

The height of the first plexiglass cylinder was 25 cm. A perforated
plexiglass plate and wire mesh were embedded on the inner wall 15 cm
away from the bottom of the particle collection device to separate the
collection and erosion simulation devices (Fig. 3). The perforated
plexiglass plate had a diameter of 19.5 cm and was filled with round
holes with a diameter of 1 cm, and the wire mesh had an aperture of 1
cm. The narrow gap between the perforated plate and the erosion
simulation device was filled with geotextiles to eliminate the boundary
effect of the particles as much as possible. Two circular holes with a
diameter of 5 cm were arranged on relative positions in the pipe wall to
form a hydraulic connection between the water tank and the erosion
simulation device during the saturation and erosion tests. The two
circular holes and the bottom of the erosion particle collection chamber
were sealed with non-woven geotextiles, and only water exchange was
allowed between the sample and the water tank (Fig. 3).

The data collection devices used in this test included the following:
pore water pressure gauge, dial indicator, and DH3816N data acquisition
and analysis system. The pore water pressure gauge was installed in each
ring of the erosion simulation device, and the gap between the conductor
of the pore water pressure gauge and the hole was sealed with waterproof
mastic to avoid the formation of a seepage channel (see Fig. 4).
Moreover, a perforated plexiglass plate with a diameter of 20 cm and
weight of 1 kg filled with round holes with a diameter of 1 cm was
placed on the upper surface of the sample. A dial indicator for
measuring the displacement was positioned at the centre of the weight to
measure the settlement displacement of the sample.

\subsubsection*{3.2 Particle Size
Distribution}

{\label{particle-size-distribution}}

The research object of this test was soil--rock mixtures with
grap-graded PSDs. The particle size at the boundary between the coarse
and fine particles was selected according to the British Standard (1990)
and the Ministry of Water Resources of the People's Republic of China
(1999). At the same time, considering the test operability and data
analysis, 2 mm was chosen as the boundary particle size of the coarse
and fine particles. Quartz sand with a particle size of 0.5--2 mm was
selected as the fine particles. Meanwhile, the maximum particle size of
the sample should not be greater than 1/5 of the test size (Ministry of
Water Resources of the People's Republic of China, 1999). Therefore,
sandstone with a particle size of 10--30 mm was selected as the coarse
material, which was obtained from a quarry near Chongqing, China. The
ratio of the coarse to fine particles was 5:60.

Based on the fine particle content, Vallejo et al. (2001) and Liu (2006)
divided the soil skeleton into coarse, coarse--fine, and fine particle
skeleton structures. In the coarse--fine particle skeleton structure,
the FC was approximately 30\%--40\% (Minh et al., 2014; Lopez et al.,
2016). Coarse particles act as soil skeletons and bear most of the load,
whereas fine particles fill the pores between the coarse particles and
transfer part of the effective stress. However, under rainfall
conditions, infiltration destroys the strong interaction between the
coarse and fine particles, and as the FC decreases, the skeleton
structure becomes unstable. Therefore, it is necessary to analyse the
changes in the skeleton structure with FC ranging from 30\% to 40\%
under the effect of infiltration. In this study, three FCs, 30\%, 35\%,
and 40\%, were selected to analyse the seepage deformation of
grap-graded soil--rock mixtures. The original particle size distribution
of the test material is shown in Fig. 5. The specific gravity of the
raw-material soil particles was regarded as approximately 2.65
g/cm\textsuperscript{3}. The height of the filled samples was 50 cm, and
the bottom area was 314.16 cm\textsuperscript{2}. Table 1 summarises the
relevant parameters of the materials used in this study.

\subsubsection*{3.3 Preparation and Test
Procedure}

{\label{preparation-and-test-procedure}}

A batch layered preparation method was used to reshape the sample
(Bendahmane et al., 2008), and a sample with uniform mixing of coarse
and fine particles was obtained. The specific steps are as follows:

\begin{enumerate}
\tightlist
\item
  The total mass of each batch of materials under dry conditions was 2.0
  kg, and the corresponding masses of the quartz sand and sandstone were
  weighed separately according to the different proportions of coarse
  and fine particles in each test.
\item
  Water with a slightly lower standard than that of the optimal water
  content was added into the quartz sand and stirred.
\item
  According to the set proportion, sandstone was added into the quartz
  and water mixture and stirred to obtain the final soil--rock mixture.
  If too many soil--rock mixtures are prepared at once, it is easy to
  cause a large separation between the fine and coarse particles at this
  stage.
\end{enumerate}

The soil--rock mixture was poured into the erosion simulation device and
stirred several times. Although some fine particles in the first batch
passed through the wire mesh into the particle collection device before
the test, there was almost no particle loss afterwards, which was also
consistent with the test phenomenon of Moffat et al. (2006). The water
inlet pipe was close to the inner wall of the water tank to ensure that
the water level slowly infiltrated the sample from bottom to top to
minimise the impact of seepage force on the stability of the fine
particles. The method of static saturation (Ministry of Water Resources
of the People's Republic of China, 1999) was adopted to ensure that the
sample remained for 12 h. Based on the observation, there were no
bubbles visible to the naked eye in the saturated sample.

This test adopted the seepage mode of the constant head (I = 0.2) and
constant load (1.0 kg) and was carried out approximately 2 h after
saturation. When the test was completed, the drain valve at the bottom
of the water tank was opened to allow all the water in the tank and
erosion simulation device to be discharged. Afterwards, different
samples were collected from the five layers of the erosion simulation
device and stored separately. After drying and weighing, the
distribution changes in each layer were analysed.

\subsubsection*{3.4 Test Results}

{\label{test-results}}

\emph{3.4.1 Pressure Gradient Variation}

{\label{pressure-gradient-variation}}

Fig. 6 shows the pressure gradient variation of the five layers in the
erosion simulation device during the test, where i\textsubscript{jk}
represents the local hydraulic gradient variation between the j-th and
k-th pore pressures.

For FC = 30\% and 35\% in Figs. 6(a) and (b), respectively,
^i\textsubscript{1-2} and ^i\textsubscript{4-5} (bottom and top of the
sample) were greater than ^i\textsubscript{2-3} and
^i\textsubscript{3-4} (middle of the sample). In particular,
^i\textsubscript{3-4} slightly changed during the test. This may be due
to the loss of fine particles and internal latent corrosion of the
first, second, fourth and fifth layers. Although the fine particles in
the third and fourth layers also migrated, they were supplemented from
the top of the sample (fifth layer). Moreover, under the coupling
mechanism of particle loss and particle supplementation, the pressure
gradient variation slightly changed. The changes in the two samples
remained stable after approximately 1500 and 2000 s, respectively. The
larger the FC, the more the particles, and the longer the stabilisation
time. Furthermore, ^i\textsubscript{1-2} and ^i\textsubscript{2-3}
(upstream of the sample) were greater than ^i\textsubscript{3-4} and
^i\textsubscript{4-5} (downstream of the sample), which is also
consistent with the test phenomenon of Jiang et al. (2018). For FC =
40\%, as shown in Fig. 6(c), ^i\textsubscript{1-2}and
^i\textsubscript{2-3} were greater than ^i\textsubscript{3-4} and
^i\textsubscript{4-5} . However, ^i\textsubscript{3-4} increased after
1200 s and continued to expand at the end of the test, which may be
because the seepage channel of the sample with FC = 40\% is more
difficult to stabilise, causing the erosion at the top of the sample to
develop to the middle part.

\emph{3.4.2 Changes in Fine Particles in Each Layer}

{\label{changes-in-fine-particles-in-each-layer}}

Fig. 7 shows the changes in fine particles in each layer after the test.
It was found that the top change of fine particles in the erosion
simulation (first layer) was the largest, followed by the bottom change
(fifth layer), and then the middle change (second, third, and fourth
layers). The greater the FC, the greater the change in fine particles
over the entire height. Moreover, the loss of fine particles in the
sample with FC = 40\% in the fifth layer reached a maximum of 5.3\%.

After fine particle migration under the action of infiltration, the
vertical section of the sample was divided into three areas according to
changes in fine particles: top, middle uniform, and bottom loss areas.
The loss of fine particles in the top loss area was the largest,
followed by the loss in the bottom loss area, and then the loss in the
middle uniform area. Because fine particles at the top of the sample
were lost to the next layer of the sample without other supplements, the
loss of fine particles was the largest. The bottom layer of the sample
was located in the downstream section of the seepage channel, and the
fine particles were prone to large loss. Meanwhile, the fine particles
in the middle of the sample were lost and received supplements from the
top of the sample; hence, the change in fine particles was the smallest.
Similar phenomena were also observed in the permeability deformation
test (Kenney et al., 1985; Marot et al., 2016; Zhong et al., 2018).

It should be mentioned that because the operation error during the test
resulted in the successful recording of the dial indicator data of only
one of the three samples, there was no comparative analysis of the top
settlement of different samples.

\subsection*{PARTICLE-SCALE NUMERICAL
TEST}

{\label{particle-scale-numerical-test}}

\subsubsection*{4.1 PFC Model}

{\label{pfc-model}}

The numerical test was simplified based on the physical test and the
calculation workload. Only ideal discontinuous grap-graded soil--rock
mixtures with fine (particle size range of 1--2 mm) and coarse particles
(particle size range of 10 mm) were considered. The particle shape used
in the simulation was spherical, and its influence was not considered.
Three types of PSD were selected, and the mass fractions of fine
particles in the three gradations were 30\%, 35\%, and 40\%. During
sample preparation and testing, the gravitational acceleration was set
to 9.81 m/s\textsuperscript{2}, and the direction was vertical and
downward.

The numerical test piece was a cylinder with a diameter of 50 mm and
height of 125 mm. The contact model in the PFC adopts a rolling
resistance linear model, and Table 2 lists the contact parameters of the
PFC particles. In Table 2, a magnitude of 1.5e8 Pa was selected for the
Young's modulus in this study. Although the elastic models of quartz and
sandstone are, in reality, greater than 1e9 Pa, Chand et al. (2012)
found that a smaller magnitude of Young's modulus does not significantly
affect the physical response of the particulate system and also allows
for a larger time step, which can considerably save computational
resources.

Presently, there are several methods of preparing accumulation particle
samples using DEM, such as the falling method (\selectlanguage{ngerman}ØREN et al., 2002), pure
geometric method (Bagi et al., 2005), radius amplification method
(Itasca Consulting Group I, 2004), and layered under compaction method
(Jiang et al., 2003). The object of this study belongs to grap-graded
PSDs under compaction, that is, natural accumulation. To simulate this
stacking state, the sample was prepared based on the works of Huang et
al. (2015) and Huang et al. (2020). The specific steps for generating
the particle accumulation are as follows:

Stage 1: Initial particles were generated. According to the given soil
particle distribution curve, particles were randomly generated in the
model by the Monte Carlo method in the wall, which is twice the height
of the target model as the initial state (Fig. 8(a)). The initial
porosity of the generated particles was set to 0.7, and it was used to
control only the number of initial particles.

Stage 2: Stable particle accumulation was obtained after gravitational
deposition. By applying gravitational field to the particles, the
particles were stably deposited under the action of self-weight, and
their porosities were measured. A servo pressure of 100 Pa was applied
on the deposited particle to compact to pressurise the model to the
preset sample height. Similar to the physical test, the sample was
divided into five layers based on height (Fig. 8(b)).

Stage 3: The sample was obtained. The lower wall was removed, and a
perforated diaphragm with a circular hole diameter of 6 mm and spacing
of 10 mm was placed (Fig. 8(c)). In this stage, the porosities of the
samples with FC = 30\%, 35\%, and 40\% were 0.39, 0.36, and 0.34,
respectively. As shown in Fig. 8(d), the sample has a coarse--fine
particle skeleton structure, and the coarse and fine particles were
stressed together.

\subsubsection*{4.2 ABAQUS model}

{\label{abaqus-model}}

ABAQUS/Standard calculation program was introduced to calculate the
seepage in the sample. The shape and size of the ABAQUS calculation grid
were consistent with the DEM samples, and there were 120 elements in
total, as shown in Fig. 9. As in the physical test, the upper boundary
condition of the model was set as the pressure inlet, the lower boundary
was set as the pressure outlet, and the surrounding boundary conditions
were set as the impermeable wall boundary to simulate one-way seepage.
The pressure gradient was set to 0.2. The fluid parameters used in this
study are listed in Table 3. All particle information (including
position, velocity, and drag force) and all contact forces and energies
within the samples were recorded every 0.1 s during each suffusion
simulation, with a total duration of 20 s. Each coupling interval was
approximately 20 min, and with a total of 200-time intervals, it took
approximately 2--3 days on a sample. The calculation was conducted on a
computer with Intel I9-9900K and 32 GB RAM.

\subsection*{NUMERICAL RESULTS}

{\label{numerical-results}}

\subsubsection*{5.1 Pressure Gradient
Variation}

{\label{pressure-gradient-variation-1}}

Fig. 10 shows the pressure gradient variations of the different samples
during seepage. The fluid element was evenly divided into five parts
according to the physical test. In the numerical test of three samples,
the pressure gradient of the specimen at different parts of the seepage
path was different. This is similar to the phenomenon and law observed
in the physical test, where \selectlanguage{english}^i\textsubscript{1-2} and
\selectlanguage{english}^i\textsubscript{4-5} (bottom and top of the sample) was greater than
\selectlanguage{english}^i\textsubscript{2-3} and \selectlanguage{english}^i\textsubscript{3-4} (middle of the sample)
and the pressure gradients of the upper part (\selectlanguage{english}^i\textsubscript{3-4} and
\selectlanguage{english}^i\textsubscript{4-5} ) was less than those of the lower part
(\selectlanguage{english}^i\textsubscript{1-2} and \selectlanguage{english}^i\textsubscript{2-3} ). For the sample with
FC = 30\%, the phenomena and laws were obvious. In addition, for the
sample with FC = 40\%, the change in pressure gradient
\selectlanguage{english}^i\textsubscript{3-4} was also obvious, which is similar to the
phenomenon observed in the physical test (Fig. 6(c)). This indicates
that the particle structure changes noticeably and is caused by the
migration of a large number of fine particles.

\subsubsection*{5.2 Transportation of Fine
Particles}

{\label{transportation-of-fine-particles}}

\emph{5.2.1 Macroresponses of the Gap-graded Samples}

{\label{macroresponses-of-the-gap-graded-samples}}

Fig. 11 shows the variation in the roof settlement and fine particle
loss during the test. Before approximately 3 s, the number of particle
losses was inversely proportional to the fine particle content. The loss
of fine particles in the sample with FC = 30\% was the largest, and the
roof settlement of the experimental sample with FC = 40\% was the
largest. However, the roof settlement of the sample with FC = 30\% was
greater than that of the sample with FC = 35\%. This may be because the
skeleton of the sample with FC = 30\% was unstable, the loss of fine
particles in the first layer was significant, and the roof settlement
was also large.

For the sample with FC = 30\%, the roof settlement and particle loss
stabilized at approximately 2.5 and 2.7 s, respectively. Meanwhile, for
the sample with FC = 35\%, the roof settlement and particle loss
stabilized at approximately 3.8 and 7.5 s, respectively. The stability
of the roof settlement and particle loss indicated that the load-bearing
skeleton was stable. For the sample with FC = 40\%, there was no
stabilisation time, and the rate of change of sedimentation and particle
loss reduced in 4 s.

\emph{5.2.2 Particle Migration}

{\label{particle-migration}}

Based on the analysis of the seepage failure process of each test, some
specific time points were selected for analysis (Fig. 11(b)). Fig. 11
depicts the internal migration of pre-stratified particles at different
times, where the black arrow represents the main infiltration channel.
The figure shows that the coarse particles moved slightly.

Before the test, the number of particles lost in the first layer was the
highest, as shown in Fig. 12(a). As the test progressed, particles in
the second and third layers entered the first layer. Afterwards,
particles in the fourth layer penetrated the third layer but did not
penetrate the second layer. Under the action of the fluid drag force,
almost all the particles in the fifth layer penetrated the fourth layer.
Fig. 12(b) is similar to Fig. 12(a). However, because of the finer
content of the sample, there were more main infiltration channels than
the samples with FC = 30\% and 35\%. After 7.5 s, the particles in the
second to fifth layers penetrated the first layer.

To further analyse the loss of fine particles in each layer, the number
of fine particles at different heights in the test process was counted.
As shown in Fig. 13(a), before approximately 3 s, the particle loss of
each layer tended to be stable except for the fifth layer, and at
approximately 17 s, the particles in the fourth layer began to enter the
third layer. The particle accumulation in the third layer caused the
number of particles to increase. Thus, Fig. 13(b) is similar to Fig.
13(a). However, no obvious particle accumulation was observed. In Fig.
13(c), all the particles in the fifth layer were lost, and the particle
loss in the other layers was not stable.

Fig. 14 shows the changes in fine particles in each layer after the
test. The vertical section of the sample can also be divided into three
areas according to changes in the fine particle content: top loss,
middle uniform, and bottom loss areas. This phenomenon is similar to
that of the physical experiments. By comparing Figs. 13(a)--(c), it was
found that fine particles with FC = 30\% emerged and accumulated in the
middle uniform area (second--fourth layers), causing the number of
particles to increase. Meanwhile, the number of particles in the middle
uniform area with FC = 35\% slightly changed, whereas the number of
particles with FC = 40\% gradually decreased. This shows that under
natural accumulation, the fine particles are 35\% and fill the pores of
large particles.

\subsubsection*{5.3 Coarse-Fine Particle Skeleton
Structure}

{\label{coarse-fine-particle-skeleton-structure}}

\emph{5.3.1 Particle-Particle Contacts}

{\label{particle-particle-contacts}}

There are five types of particle contact: fine--fine, coarse--coarse,
coarse--fine, coarse--wall, and fine--wall contacts. Owing to the loss
of fine particles, the fine--fine and coarse--coarse contacts of
grap-graded soil--rock mixtures change most significantly during seepage
failure. Fig. 15 shows the changes in the fine--fine and coarse--coarse
contacts.

From Fig. 15, the number of fine--fine and coarse--coarse contacts of
the sample with FC = 30\% remained stable during the test, indicating
that although there was a certain loss of fine particles under this fine
particle content, it did not affect the change in the load-bearing
skeleton. For the sample with FC = 35\%, before approximately 5 s, the
number of fine--fine contacts decreased with the loss of fine particles,
whereas the number of coarse--coarse contacts increased, showing that
the coarse particles in the load-bearing skeleton gradually bear the
load of the lost fine particles. After approximately 5 s, the number of
contacts did not change. For the sample with FC = 40\%, the number of
fine--fine contacts decreased throughout the test, whereas the number of
coarse--coarse contacts increased. At the end of the test, the number of
contacts did not remain stable, resulting in the continuous loss of fine
particles and the number of fine--fine contacts being less than that of
samples with FC = 30\%.

Fig. 16 shows the axial force nephogram of fine--fine and coarse--coarse
contacts in the load-bearing skeleton after the test. In Fig. 16(a), the
area in which the fine particles accumulate between the pores of the
coarse particles is marked in red. It is obvious from the figure that
the number of fine--fine contacts in the sample with FC = 35\% was the
largest, and the stress distribution was relatively large. In Fig.
16(b), for the sample with FC = 30\%, the force chain of the
coarse--coarse contact was more complete from the top to bottom of the
sample. For the sample with FC = 35\%, the force chain of the
coarse--coarse contact was missing in the middle from the top to bottom
of the sample, that is, the dense area of fine particles in Fig. 16(a).

\emph{5.3.2 Strong Force Chain Buckling}

{\label{strong-force-chain-buckling}}

Because of the loss of fine particles, the original load-bearing
skeleton collapsed, and a new skeleton eventually formed. The collapse
of the original load-bearing skeleton was mainly induced by strong force
chain buckling. To further understand the evolution characteristics of
the packing of particles of different load-bearing skeletons, the upper
and lower parts (the junction of the first and second layers and the
junction of the third and fourth layers) in each test were selected, as
shown in Fig. 17. The strong force chains and chained particles are
marked in red. The contact forces between the particles were larger than
the average value. Fig. 18 shows the evolution characteristics of the
force chain for samples with different FCs at the selected strong force
chains. For ease of display, the particles were reduced once.

In Fig. 18(a), the axial force decreased in the initial 1 and 2 s of the
samples with FC = 30\% and FC = 35\%, respectively. Afterwards, the fine
particles in the fifth layer gradually flowed downwards and became lost.
Moreover, the interaction of coarse particles between the fourth and
fifth layers became increasingly obvious, and the strength of the force
chain continuously increased. The samples with FC = 30\% and 35\%
remained stable (approximately 16 and 13 N) after 15 s. In the sample
with FC = 40\%, there is no force chain between the selected coarse
particles in the initial load-bearing structure (the axial force was 0).
With the development of seepage failure and the loss of fine particles,
the axial force of the force chain between the selected coarse particles
increased. After 13 s, the axial force remained relatively stable, that
is, 16--18 N.

In Fig. 18(b), the axial force decreased in the initial 2 and 4 s of the
samples with FC = 30\% and 35\%, respectively. Subsequently, the axial
force remained stable. The time was the same as the stable time of
particle loss in Fig. 11(b), indicating that the loss of fine particles
stopped and the load-bearing structure of the skeleton remained stable.
In the sample with FC = 40\%, the evolution was the same as that in Fig.
18(a). After 13 s, the axial force remained relatively stable, that is,
26--28 N.

\subsection*{DISCUSSION}

{\label{discussion}}

The relative content of fine particles is the key property affecting
stress transfer through grap-graded granular materials. The selected FC
in this test was 30\%--40\%, which is a coarse--fine particle skeleton
structure. The sample with FC = 30\% tended to have a coarse particle
skeleton structure (fine particles underfilled the voids between the
coarse particles), whereas the sample with FC = 40\% tended to have a
fine particle skeleton structure (fine particles overfilled the voids
between the coarse particles).

Under the action of vertical infiltration and gravity, the change in
pressure gradient in the middle area was less than that in the upper and
lower areas (Fig. 10). However, as the FC increased, the pressure
gradient in the middle area changed considerably. Because fine particles
play an important role in the load-bearing structure, if the fine
particles are dragged and migrated, the structure will change
significantly, causing the pressure gradient to change with time.

Fine particles were lost at the top and bottom of the three samples;
however, as the FC increased, the particles in the middle soil changed
from a certain accumulation to stability to particle loss (Figs. 12 and
13). Fine particles reasonably fill the voids between coarse particles
and can maintain better stability under the action of the same
infiltration force. Obviously, the finer the particles, the greater the
force shared by the fine particles. Therefore, in the process of fine
particle loss, the samples with FC = 30\% and 35\% were relatively more
stable during the test process. For the sample with FC = 35\%, fine
particles bore more load contribution (with a finer particle contact
number, as shown in Fig. 16(a)). For the samples with FC = 40\%, the
initial contact force of the strong link was 0 (Fig. 18), which
indicates that there was poor contact between the coarse particles such
that the coarse particles became an ``island'' and only a part of the
fine particles around the coarse particles bore the load. Meanwhile, the
coarse particles bore more load with the loss of fine particles.

By comparing the phenomena and results of the physical and numerical
tests, there is a high similarity in the variation process of the
pressure gradient and fine particles in different layers. This shows
that the numerical model can better reproduce meso phenomena such as
particle migration, permeability evolution, and particle recombination
in the physical test. The coupling model of PFC3D and ABAQUS can be used
as an effective means of studying the mechanism of seepage deformation.
However, owing to the defect in calculation efficiency of the numerical
software, the method needs to be further improved to reflect the
irregular characteristics of particles and the influence of micro
disturbances in real situations through GPU/CPU parallel technology.

\subsection*{CONCLUSIONS}

{\label{conclusions}}

In this study, a permeability test for the stability of gap-graded
soil--rock mixtures considering particle loss was conducted, and the
physical process of seepage deformation of soil--rock mixtures with
different FCs was numerically simulated using the coupling model of
PFC3D and ABAQUS. The results showed that this coupling model has a good
effect on the simulation of the particle-scale movement under seepage.
Particle migration causes a change in the soil particle structure and
also induces seepage deformation. Hence, the key concluding remarks are
as follows:

\begin{enumerate}
\tightlist
\item
  During the seepage process, fine particle loss occurred at different
  sample heights. After seepage, the spatial distribution of fine
  particle loss along the height direction can be divided into three
  areas: top, middle uniform, and bottom loss areas. The pressure
  gradient in the top loss area was smaller than that in the bottom loss
  area during the test.
\item
  The ``island'' effect of coarse particles, which is caused by
  excessive fine content and makes the fine particles bear more load,
  was eliminated with the loss of fine particles. The fine particle
  content remained relatively stable under the actions of infiltration
  and vertical gravity. In this preset working condition of coarse and
  fine particle diameters, setting FC to 35\% may be the best way to
  fill the voids between the coarse particles.
\item
  The change processes of the pressure gradient, transportation of fine
  particles, and the load-bearing structure were recorded, which clearly
  and intuitively expresses the migration trend of particles in the
  seepage process. Moreover, the numerical simulation results were
  highly consistent with the physical test results. This shows that the
  coupling model of PFC3D and ABAQUS can capture both micro- and
  macro-characteristics of particles and fluids and has practical
  significance in simulating the particle-scale seepage deformation
  process of grap-graded soil--rock mixtures.
\end{enumerate}

\subsection*{DECLARATION OF CONFLICTING
INTERESTS}

{\label{declaration-of-conflicting-interests}}

The author(s) declared no potential conflicts of interest with respect
to the research, authorship, and/or publication of this article.

\subsection*{ACKNOWLEDGMENTS}

{\label{acknowledgments}}

The author(s) disclosed receipt of the following financial support for
the research, authorship, and/or publication of this article: This paper
gets its funding from projects (Grant No. CYB21031) supported by the
graduate research and innovation foundation of Chongqing, China;
(LNTCCMA-20200103; LNTCCMA-20210107) supported by the Key Laboratory of
New Technology for Construction of Cities in Mountain Area. The authors
wish to acknowledge these supports.

\subsection*{REFERENCES}

{\label{references}}

Bagi K (2005) An algorithm to generate random dense arrangements for
discrete element simulations of granular assemblies. \emph{Granular
Matter} 7(1), 31-43.

Bendahmane F, Marot D and Alexis A (2008) Experimental parametric study
of suffusion and backward erosion. \emph{Journal of geotechnical and
geoenvironmental engineering} 134(1), 57-67.

BSI (1990). BS1377:1990: Methods of test for soils for civil engineering
purposes, Part 1: General requirements and sample preparation.

\emph{British: Authority of the Board of BSI Press.}

Carman PC (1939) Permeability of saturated sands, soils and
clays.\emph{The Journal of Agricultural Science} 29(2), 263-273.

Chand R, Khaskheli MA, Qadir A, et al. (2012) Discrete particle
simulation of radial segregation in horizontally rotating drum: Effects
of drum-length and non-rotating end-plates. \emph{Physica A-Statistical
Mechanics and Its Applications} 391(20), 4590-4596.

Cheng K, Wang Y and Yang Q (2018) A semi-resolved CFD-DEM model for
seepage-induced fine particle migration in gap-graded
soils.\emph{Computers and Geotechnics} 100, 30-51.

Cheng K, Zhang C, Peng K, et al. (2021) Un-resolved CFD-DEM method: an
insight into its limitations in the modelling of suffusion in gap-graded
soils. \emph{Powder Technology} 381: 520-538.

Cui Y, Nouri A, Chan D, Rahmati E (2016) A new approach to DEM
simulation of sand production. \emph{J Petrol Sci Eng} 147, 56-67.

Difelice, R (1994) The voidage function for fluid-particle interaction
systems. \emph{Int. J. Multiphase Flow} 20(1), 153-159.

Fan SY, Song ZP, Zhang YW, et al. (2020) Case Study of the Effect of
Rainfall Infiltration on a Tunnel Underlying the Roadbed Slope with Weak
Inter-Layer. \emph{KSCE Journal of Civil Engineering} 24(5), 1607-1619.

Goodarzi M Kwok CY and Tham LG (2014) A continuum-discrete model using
Dacry's Law: formulation and verification. \emph{Int J Numer Anal Meth
Geomech} 29(3), 327-342.

Grosta G and Prisco C (1999) On slope instability induced by seepage
erosion. \emph{Canadian Geotechnical Journal} 36 (6), 1056-1073.

Guo Y and Yu XB (2017) Comparison of the implementation of three common
types of coupled CFD-DEM model for simulating soil surface
erosion.\emph{International Journal of Multiphase Flow} 91, 89-100.

Huang QF, Zhan ML, Sheng JC, et al. (2015). Numerical method to generate
granular assembly with any desired relative density based on
DEM.\emph{Chinese Journal of Geotechnical Engineering} 37(03), 537-543.

Huang Y, Luo YL and Zhang XG. (2020). Numerical method to generate
particle packages of internally unstable soils with arbitrary relative
density. \emph{Water resources and power} 38(04), 126-129.

Itasca Consulting Group I. (2004) User's manual: PFC3D-particle flow
code in 3 dimensions. \emph{Minneapolis: Itasca Consulting Group} .

Itasca Consulting Group Inc, 2018. PFC3D 5.00.
Documentation.\emph{Minneapolis: Itasca Consulting Group Inc} .

Jiang MJ, Konrad JM and Leroueil S (2003) An efficient technique for
generating homogeneous specimens for DEM studies. \emph{Computers and
geotechnics} 30(7), 579-597.

Jiang ZM, Yuan T, Liu DQ, et al. (2018). Mesoscopic numerical tests on
seepage failure characteristics of coarse grained soils. \emph{Chinese
Journal of Geotechnical Engineering} 40(4), 752-758. (in Chinese)

Johan V and Zou QP (2013) Investigation of slope instability induced by
seepage and erosion by a particle method. \emph{Computers and
Geotechnics} 48: 9-20.

Kenney TC and Lau D (1985) Internal stability of granular
filters.\emph{Can Geotech J} 22(2), 215-225.

Kozeny J (1927) Ueber kapillare Leitung des Wassers im
Boden.\emph{Sitzungsber Akad. Wiss, Wien} 136(2a), 271-306.

Lei XQ, Yang ZJ, He S, et al. (2017) Numerical investigation of
rainfall-induced fines migration and its influences on slope
stability.\emph{Acta Geotech} 12, 1431-1446.

Liu J (2006) Seepage control of earth-rock dams theoretical
basis.\emph{Beijing: Chian Water \& Power Press} . (in Chinese)

Liu YJ, Wang LZ, Hong Y, et al. (2020) A coupled CFD-DEM investigation
of suffusion of gap graded soil: Coupling effect of confining pressure
and fines content. \emph{International Journal for Numerical and
Analytical Methods in Geomechanics} 44(18), 2473-2500.

Lopez R, De F, Silfwerbrand J, et al. (2016) Force transmission and soil
fabric of binary granular mixturese. \emph{Geotechnique} 66(7), 578-583.

Marot DA, Rochim HH, et al. (2016) Assessing the susceptibility of gap
graded soils to internal erosion:proposition of a new experimental
methodology. \emph{Natural. Hazards} 83(1), 365-388.

Minh NH, Cheng YP and Thornton C (2014) Strong force networks in
granular mixtures. \emph{Granular Matter} 16(1), 69-78.

Ministry of Water Resources of the People's Republic of China (1999).
GB/T 50123-1999. \emph{China Planning Press} .

Mitchell JK and Soga K (2005) Fundamentals of Soil Behavior. \emph{New
York: John Wiley \& Sons} , 2005.

Moffat R, Ricardo A, Fannin, et al. (2006) A large permeameter for study
of internal stability in cohesionless soils. \emph{Geotechnical Testing
Journal} 29(4), 273-279.

Oren P, Bakke S (2002) Process based reconstruction of sandstones and
prediction of transport properties. \emph{Transport in Porous
Media}46(2/3), 311-343.

Ouyang M and Takahashi A (2015) Influence of initial fines content on
fabric of soils subjected to internal erosion. \emph{Canadian
Geotechnical Journal} 53(2), 1360-1361.

Peters B, Baniasadi M, Baniasadi M, et al. (2019) XDEM multi-physics and
multi-scale simulation technology: Review of DEM-CFD coupling,
methodology and engineering applications -
ScienceDirect.\emph{Particuology} 44, 176-193.

Shi DM, Wang WL, Jiang, GY, et al. (2016) Effects of disturbed landforms
on the soil water retention function during urbanization process in the
Three Gorges Reservoir Region, China. \emph{Catena} 144: 84-93.

Skempton AW and Brogan JM (1994) Experiments on piping in sandy
gravels.\emph{Geotechnique} 44(3), 449-460.

Tang Y, Chan DH, Zhu DZ (2017) A coupled discrete element model for the
simulation of soil and water flow through an orifice.\emph{International
Journal for Numerical and Analytical Methods in Geomechanics} 41(14),
1477-1493.

Thevanayagam S, Shenthan T, Mohan S, et al. (2002) Undrained fragility
of clean sands, silty sands, and sandy silts. \emph{Journal of
Geotechnical and Geoenvironmental Engineering} 128(10), 849-859.

Tian DL, Xie Q, Fu X, et al. (2020) Experimental study on the effect of
fine contents on internal erosion in natural soil
deposits.\emph{Bulletin of Engineering Geology and the Environment} 79,
4135-4150.

Tsuji Y, Kawaguchi T and Tanaka T (1993) Discrete particle simulation of
two-dimensional fluidized bed. \emph{Powder Technol} 77, 79-87.

Vallejo LE (2001) Interpretation of the limits in shear strength in
binary granular mixtures. \emph{Canadian Geotechnical Journal} 38(5),
1097-1104.

Wang Y, Chai JR, Xu ZG, et al. (2020) Numerical Simulation of the
Fluid-Solid Coupling Mechanism of Internal Erosion in Granular
Soil.\emph{Water} 12, 137.

Xiong H, Yin ZY, Zhao JD dong, et al. (2021) Investigating the effect of
flow direction on suffusion and its impacts on gap-graded granular
soils. \emph{Acta Geotechnica} 16, 399-419.

Xu BH and Yu AB (1997) Numerical simulation of the gas-solid flow in a
fluidized bed by combining discrete particle method with computational
fluid dynamics. \emph{Chemical Engineering Science} 52(16), 2785-2809.

Xu WJ, Hu RL and Tan RJ (2007) Some geomechanical properties of
soil-rock mixtures in the Hutiao Gorge area, China.
\emph{Geotechnique}57(3), 255-264.

Yang Q, Cheng K, Wang, Y, et al. (2019) Improvement of semi-resolved
cfd-dem model for seepage-induced fine-particle migration: eliminate
limitation on mesh refinement. \emph{Computers and Geotechnics} 110,
1-18.

Yin Y, Cui Y, Tang Y, et al. (2021) Solid-fluid sequentially coupled
simulation of internal erosion of soils due to seepage. \emph{Granular
Matter} 23(2).

Zhang LL, Zhang J and Zhang LM (2011) Stability analysis of
rainfall-induced slope failure: A review. \emph{Geotechnical
Engineering} 164(5), 299-316.

Zhang S, Tang HM, Zhan, HB, et al. (2015) Investigation of scale effect
of numerical unconfined compression strengths of virtual
colluvial-deluvial soil-rock mixture. \emph{International Journal of
Rock Mechanics and Mining Sciences} 77, 208-219.

Zhong CH, LE VT, Bendahamane F, et al. (2018) Investigation of Spatial
Scale Effects on Suffusion Susceptibility. \emph{Journal of Geotechnical
and Geo-environmental Engineering} 144(9), 04018067.

Zou YH, Chen C, Zhang LM, et al. (2020) Simulating Progression of
Internal Erosion in Gap-Graded Sandy Gravels Using Coupled
CFD-DEM.\emph{International Journal of Geomechanics} 20(1), 04019135.

Zou YH, Chen Q, Chen XQ, et al. (2013) Discrete numerical modeling of
particle transport in granular filters. \emph{Computers and
Geotechnics}47, 48-56.

Table 1. Physical properties of the specimen\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
\begin{minipage}[t]{0.16\columnwidth}\raggedright\strut
Number of samples\strut
\end{minipage} & \begin{minipage}[t]{0.16\columnwidth}\raggedright\strut
\emph{d}\textsubscript{10} (mm)\strut
\end{minipage} & \begin{minipage}[t]{0.16\columnwidth}\raggedright\strut
\emph{d}\textsubscript{30} (mm)\strut
\end{minipage} & \begin{minipage}[t]{0.16\columnwidth}\raggedright\strut
\emph{d}\textsubscript{60} (mm)\strut
\end{minipage} & \begin{minipage}[t]{0.16\columnwidth}\raggedright\strut
\emph{C}\textsubscript{u}(d\textsubscript{60}/d\textsubscript{10})\strut
\end{minipage} & \begin{minipage}[t]{0.16\columnwidth}\raggedright\strut
\emph{C}\textsubscript{c}()\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
FC = 30\%\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
0.82\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
3\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
16.12\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
19.75\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
0.68\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
FC = 35\%\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
0.78\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
1.65\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
15.49\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
19.87\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
0.23\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
FC = 40\%\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
0.74\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
1.28\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
14.76\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
19.99\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
0.15\strut
\end{minipage}\tabularnewline
\bottomrule
\end{longtable}

Table 2. Calculation parameters of the PFC3D model\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
\textbf{Parameter} & \textbf{Value} & \textbf{Value}\tabularnewline
& \textbf{Particle} & \textbf{Wall}\tabularnewline
Particle density (kg/m\textsuperscript{3}) & 2630 & -\tabularnewline
Effective modulus (Pa) & 1.5e8 & 1.0e9\tabularnewline
Normal-to-shear stiffness ratio & 1.0 & 1.0\tabularnewline
Contact friction coefficient & 0.5 & 0.01\tabularnewline
Normal critical damping ratio & 0.2 & 0.2\tabularnewline
Shear critical damping ratio & 0.2 & 0.2\tabularnewline
Rolling resistance coefficient & 0.1 & 0.1\tabularnewline
Time step & 1e-7 & -\tabularnewline
\bottomrule
\end{longtable}

Note: Based on the existing test results, Mitchell et al. (2005)
suggested that the friction coefficient of quartz sand is 0.5. The
common anti-rotation coefficient of the particle is taken as 0.1, based
on the literature (Liu et al., 2020).

Table 3. Calculation parameters during the coupling process (Note: the
porosity, particle applied force, and drag force are updated every 20
cycles in each coupling interval.)\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
\textbf{Parameter} & \textbf{Value}\tabularnewline
Fluid density (kg/m\textsuperscript{3}) & 1e3\tabularnewline
Fluid dynamic viscosity (Pa\selectlanguage{ngerman}·s) & 1e-3\tabularnewline
Gravitational acceleration (m/s) & 9.8\tabularnewline
Pressure gradient difference & 0.2\tabularnewline
DEM timestep (s) & 1e-7\tabularnewline
Coupling interval (s) & 0.1\tabularnewline
CFD update interval & 20\tabularnewline
Total simulation time (s) & 200\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-1-Calculation-sequence-for-the-fluid-solid-interaction-process/Figure-1-Calculation-sequence-for-the-fluid-solid-interaction-process}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-2-Sketch-of-the-coarse-grid-approach/Figure-2-Sketch-of-the-coarse-grid-approach}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-3-Representation-of-sample-saturation-device,-erosion-simulation-device-with-constant-stead,-and-particle-collection-device/Figure-3-Representation-of-sample-saturation-device,-erosion-simulation-device-with-constant-stead,-and-particle-collection-device}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-4-Layout-of-data-collection-devices/Figure-4-Layout-of-data-collection-devices}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-5-Particle-size-distribution-of-the-test-material/Figure-5-Particle-size-distribution-of-the-test-material}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-6-Pressure-gradient-variation-of-samples-(a)FC=30-;-(b)FC=35-;-(c)FC=40-/Figure-6-Pressure-gradient-variation-of-samples-(a)FC=30-;-(b)FC=35-;-(c)FC=40-}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-7-Change-of-fine-particles-(note-it-is-assumed-that-the-total-amount-of-fine-particles-before-the-test-is-m-and-after-the-test-is-mj,-so-the-change-of-fine-particles/Figure-7-Change-of-fine-particles-(note-it-is-assumed-that-the-total-amount-of-fine-particles-before-the-test-is-m-and-after-the-test-is-mj,-so-the-change-of-fine-particles}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-8-Model-setup-of-the-sample-with-FC=-35--(a)-particles-in-the-stage-1;-(b)-particles-in-the-stage-2;-(c)-particles-in-the-stage-3;-(d)-contacts-in-the-the-stage-3/Figure-8-Model-setup-of-the-sample-with-FC=-35--(a)-particles-in-the-stage-1;-(b)-particles-in-the-stage-2;-(c)-particles-in-the-stage-3;-(d)-contacts-in-the-the-stage-3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-9-Fluid-computing-element/Figure-9-Fluid-computing-element}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-10-Pressure-gradient-variation-of-numerical-results-(a)-FC=30-;-(b)-FC=35-;-(c)-FC=40-/Figure-10-Pressure-gradient-variation-of-numerical-results-(a)-FC=30-;-(b)-FC=35-;-(c)-FC=40-}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-11-Macroresponses-of-the-gap-graded-samples-of-different-FC-(a)Roof-subsidence-evolution;-(b)-Evolution-of-fine-particle-loss/Figure-11-Macroresponses-of-the-gap-graded-samples-of-different-FC-(a)Roof-subsidence-evolution;-(b)-Evolution-of-fine-particle-loss}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-13-Number-of-fine-particles-at-different-heights-(a)-FC=30-;-(b)-FC=35-;-(c)-FC=40-/Figure-13-Number-of-fine-particles-at-different-heights-(a)-FC=30-;-(b)-FC=35-;-(c)-FC=40-}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-14-Vertical-profile-of-mass-loss-of-fine-particles/Figure-14-Vertical-profile-of-mass-loss-of-fine-particles}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-15-Change-of-contacts-(a)Fine-fine-contact;-(b)Coarse-coarse-contact/Figure-15-Change-of-contacts-(a)Fine-fine-contact;-(b)Coarse-coarse-contact}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-16-Load-bearing-skeleton-after-test-(a)Fine-fine-contact;-(b)Coarse-Coarse-contact/Figure-16-Load-bearing-skeleton-after-test-(a)Fine-fine-contact;-(b)Coarse-Coarse-contact}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-17-Particles-profile-and-contact-force-chain-of-local-packing-after-test-(a)-FC=30-;-(b)-FC=35-;-(c)-FC=40-/Figure-17-Particles-profile-and-contact-force-chain-of-local-packing-after-test-(a)-FC=30-;-(b)-FC=35-;-(c)-FC=40-}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-18-Evolution-characteristics-of-the-force-chain-for-the-samples-with-different-FC-(a)-upper-part;-(b)lower-part/Figure-18-Evolution-characteristics-of-the-force-chain-for-the-samples-with-different-FC-(a)-upper-part;-(b)lower-part}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

