# Variation Characteristics of Mass-Loss Rate in Dynamic Seepage System of the Broken Rocks

Authors: Luzhen Wang; Hailing Kong
Journal: Geofluids
 (2018)

---

## Abstract

When the collapse column and its adjacent rocks in complex geological structures are disturbed by mining, concomitant fine particle migration, mass loss, and porous structure variation during the water seepage process in broken rocks are the inherent causes for collapse column activation and water inrush. Studying the time-varying characteristics of the mass-loss rate in the dynamic seepage system of the broken rocks is of theoretical importance for the prevention of water inrush from the collapse columns. In this study, the seepage tests of the broken mudstone were carried out using the patented pumping station seepage method, and the time-varying function of the mass-loss rate was generalized. Then, the optimal coefficients in the function of mass-loss rate were computed using the genetic algorithm. At last, the mass-loss rate in the dynamic seepage system of the broken rocks with consideration of the acceleration factor was calculated using Lagrange discrete element method. The results showed that (1) the mass-loss rate was expressed as a time-dependent, exponential function with its coefficient related to the porosity, and its time-varying characteristics were affected by gradation; (2) the time-varying curves with Talbol power exponents less than 0.6 had a rapid change stage and a slow change stage, while the time-varying curves with Talbol power exponents greater than 0.6 had an initial gradual change stage, a rapid change stage and a slow change stage; (3) at the early seepage stage, the mass-loss rate decreased with Talbol power exponent increasing; and (4) after long time seepage, the mass-loss rate was close to zero and unrelated to Talbol power exponent, and the porous structure in broken rocks remained stable with its porosity close to a certain stable value.


---

## Body

## 1. Introduction

At present, China is still a country with coal as its major energy source. With the gradual exhaustion of its shallow resources, deep excavation of coal resources has become the norm with the mining depth of individual coal mines of down to 1500 m [1]. However, deep coal seams are characteristic of a very complex geological structure, such as faults and collapse columns. Thus, their mining is often accompanied with severe water inrush, bringing serious threats to the safe production of coal mines [2, 3].As a good underground water passageway, the collapse column is one of the main culprits for frequent water inrush in deep mine exploitation. Taking widely developed collapse columns in the North China coalfield as an example, they are broken rocks with variable porous and fracture structures formed by mixing various broken rock fragments of different sizes, cementing together through muddy and siliceous thin layers on their surfaces and compacting under the self-gravity of collapse columns or ground pressure [4]. Many collapse columns with their diameters of up to tens of meters often run through multiple strata. Due to mining disturbance, internal pores and fractures of rocks develop and expand, leading to changes in the stress field and the seepage field inside themselves and their adjacent rocks, subsequently, evident alteration of the network of cracks in the confined rocks and the porous structure of the broken rocks, and finally, the formation of water seepage passages.Acted by groundwater pressure, their interior broken rocks undergo dissolution, erosion, and abrasion, and the resultant fine particles migrate and disappear under the entrainment of water seepage. When the migration and disappearance develops to a certain degree, the regional and connected pipe flow passages are formed inside the broken rocks and the confined underwater penetrates the roof and floor of coal seams through these passages, resulting in the occurrence of water inrush disasters.Many related studies [4–15] have shown that during water inrush from the collapse column, coal mining makes the roof and floor strata of the coal seams generate a large number of fractures whose expansion and connection with the collapse column leads to water inrush, that is, indirect water inrush of the broken rocks. Clearly, these studies ignored the structural features of the collapse column itself. According to the on-site observation data, a small amount of mud and sand flows out before water inrush from the collapse column, and a large amount of rock fragments and coal debris pours out during water inrush [4]. Therefore, the migration and loss of internal fine particles during the seepage in the broken rocks is the intrinsic cause for collapse column activation and water inrush. This kind of phenomenon is frequently observed in the geotechnical field.In recent years, many researchers have applied theoretical analysis, numerical simulation, and laboratory tests to study mass loss phenomena occurring in geotechnical engineering. Kenney and Lau [16] experimentally investigated the migration of mobile fine particles driven by seepage in the framework of pores and proposed the method for determining the mass loss of mobile fine particles using the particle gradating curve. Sterpi [17] established a finite difference model for these particles utilizing the continuity equation of mobile fine particles and the empirical formula of mass loss and pressure gradient of these particles. Cividini and Gioda [18] introduced the effect of seepage on the migration and mass loss of mobile fine particles and described the erosion and migration mechanism of mobile fine particles using the finite element method. Fox et al. [19] studied the erosion of groundwater on the riverbank and found that the erosion-induced mass loss of particles is related to the seepage velocity. Annamaria et al. [20] analyzed the various characteristics of the time-related fine particle density at various points in the effective seepage channel with coarse particles as the immobile framework using the finite element method. Fujisawa et al. [21] numerically simulated and reconstructed the piping phenomenon of soil grains due to erosion and framework migration and found that the amount of soil particles migrated was affected by the saturation of soil matrix and the distribution of particle sizes. Chang and Zhang [22] proposed well-graded soil and gap-graded soil to study the internal stability and found the ability of the coarse fraction of a soil to prevent the loss of its fine fraction due to seepage flow. Ke and Akihiro [23] conducted a series of seepage tests and pointed out that the hydraulic gradient would drop with the progress of suffusion and the eroded away fine particles insrease with the  initial  amount increasing. Chen et al. [24–26] through their experimental studies found that the variation of the postpiping sand outflow over time satisfies Boltzmann’s nonlinear relationship. Ma et al. [27] through the gas seepage properties of crushed coal specimens measurement found that (1) particle crushing during compaction was the main reason to increase small-size, indicating that gas seepage properties were strongly influenced by the particle size and axial displacement, and (2) the porosity decreased with the weight loss of larger particle size.All of the above studies consistently showed that coarse particles served as the immobile framework while fine particles were the movable particles and could migrate and drain away along with seepage. These studies found that the migration and loss of mobile fine particles were affected by many factors, such as seepage type [16, 18, 24–26, 28], particle size graduation [19, 21–23, 27, 29–32], seepage velocity [20], and pore water pressure [17, 32–34]. However, they did not consider the porous structure inside the framework particles, the changes concurring with the migration and dissipation of fine particles, and the effect of time-varying porosity on fine particle loss. In essence, with the migration and loss of fine particles, the porosity of the framework structure changes over time. At the same time, the framework structure with time-varying porosity also affects the migration and dissipation of its internal fine particles. Therefore, the mass-loss rate of fine particles inside the framework of coarse particles has time-dependent, nonlinear characteristics, and it is a critical physical parameter in the seepage system. Studying the varying characteristics of the mass-loss rate is the basics to research the water-inrush incident in the collapsed column.In this study, we explored the time-varying and nonlinear characteristics of the dynamic seepage system of broken rocks and examined the varying behavior of the mass-loss rate. Based on the results measured in the lab, we generalized a time-dependent, nonlinear function of the mass-loss rate and obtained the optimal coefficient values of the function using the genetic algorithm. We then introduced a time-related mass-loss rate function by reasonably selecting the acceleration factor to establish the actual collapse column model. At last, we applied the discrete element method to design the program of mass-loss rate response in the dynamic seepage system of broken rocks to compute and analyze the mass-loss rate and the spatiotemporal distribution of fine particles in the broken rocks of the actual collapse column.In order to overcome the shortcomings that previous researchers often used, the pre-set discontinuous particle size gradation, which could not be acted as the variable in numerical and theoretical calculations, we first applied Talbol continuous grading formula to grade the broken rock samples and Talbol power exponent to describe the impact of particle size distribution on the mass-loss rate in the dynamic seepage system of broken rocks. In order to avoid blindly selecting the reference value of rational coefficient due to lack of experiences, we then adopted the genetic algorithm to calculate the optimal coefficients of the mass-loss rate function.


## 2. Experiments

### 2.1. Experimental Method and Principle

The patented pumping station seepage method was employed in our experiment on the mass-loss-concurring water seepage in broken rocks. Figure1 shows the experiment apparatus and Figure 2 shows the structure of its core device: an open permeameter.Figure 1
Schematic of the experimental setup.Figure 2
Schematic and photo of the open permeameter.The permeameter filled with broken rock sample is first installed on the electronic material testing machine. Then the porosity of the rock sample,φ∗, at its natural placement state can be obtained according to the following formula:
(1)φ∗=V∗−VzV∗=1−Mπa2h∗ps,where V∗ is the volume of the broken rock sample at its natural placement state, VZ is the volume of the compacted rock sample with the same mass, M is the total mass of the broken rock sample, ρs is the mass density of the broken rock sample, a is the inner radius of the cylinder, and h∗ is the height of the rock sample at its natural placement state.Then, an axial load is applied on the broken rock sample using the electronic material testing machine so as to control the amount of compaction and measure the height of the loaded rock sample,h0. Putting h0 into the following formula can find the initial porosity before seepage:
(2)φ0=V0−VzV0=1−Mπa2h0ρs,where V0 is the volume of the rock sample after it is forced to move a certain displacement along the axial direction under loading and h0 is the corresponding height.After the rock sample is immersed in water and saturated for half an hour, the plunger pump is turned on to push water in the water tank 1 continuously flowing crossing the rock sample from the bottom of the permeameter along the pipelines. As the water continuously flows crossing the rock sample for a long time, original and secondary fine particles in the broken rock sample migrate under the entrainment of water seepage, overflowing into the tray through the overflow tank of the permeameter and passing the hose into the vibrosieve.In the experiment, the filtration method is adopted to collect lost fine particles in a timely and stepwisely manner. During the test, the lost fine particles are collected once at a regular timeΔtii=1,2,3,…,n, with a total collected times note. After drying, the particles with a mass of Δmii=1,2,3,…,n are weighed out from the lost fine particles. Obviously, the total lost mass in the rock sample during each time period is
(3)m~=∑Δmi,and the mass-loss rate within unit time and unit volume, qi, is
(4)qi=ΔmiV0Δti.In the seepage process, with the continuous migration and wastage of fine particles in the broken rock sample, the porosity of the rock sample will also continuously change over time. The current porosity of the broken rock sample at any time during the test is calculated from the lost mass as follows:(5)φi=1−M−∑Δmiπa2h0ρs,i=1,2,3,…,n,Obviously, the time-variation features of the mass-loss rateqi and the porosity φi of the broken rock sample can easily be found with the migration and loss of fine particles in the experiment. Through the pressure sensors and flow sensors installed inside the pipeline and their connected data acquisitors and computer, the pore water pressure and flow can be timely monitored and acquired.


### 2.2. Experimental Scheme

The collapse column body often contains up to 70% of clay minerals as mudstone [35]. Therefore, mudstone is taken as the test samples with measured average uniaxial compressive strength of 9.85 MPa.In order to explore the effects of both the particle sizes (gradation) and the compaction degree (initial porosity) on the time-varying characteristics of the mass-loss rate of the broken rocks in the long-term seepage, the following experimental scheme is designed.Before the experiment, all the broken mudstone particles are sieved into 8 groups with their particle sizes as 0–2.5 mm, 2.5–5 mm, 5–8 mm, 8–10 mm, 10–12 mm, 12–15 mm, 15–20 mm, and 20–25 mm and their gradations are calculated according to the following Talbol continuous grading formula:(6)Pd=MdM×100%=dDn×100%,where P is the percentage of the particles mass whose diameter is less than d, Md is the mass of broken rock sample with particle size no greater than d, M is the total mass of broken rock sample, d is the particle diameter of the current broken rock sample, D is the maximum particle size of broken rock sample, and n is Talbol power exponent.Taking into account all the human, financial, and material resources and other factors consumed in the experiment, the values of Talbol power exponent,n, for the test samples are assigned as 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0, respectively. Based on this, the results of other gradation samples can be obtained by numerical computation, so as to not only save the cost of the test but also obtain more calculated data than the tested results.In the test, the total mass of mudstone filled in each cylinder is 2000 g. Figure3 shows the mass distribution of experimental samples with 8 kinds of particle sizes in 10 gradations. Figure 4 shows the test samples mixed uniformly at n=0.4 and n=1.0 according to the Talbol continuous grading formula.Figure 3
Mass distributions of test samples with various particle sizes in 10 gradations.Figure 4
The test samples mixed according to Talbol continuous grading formula.
(a)
n=0.4
(b)
n=1.0In the experiment, the broken rock samples are compacted to 5 different heights, as shown as the compaction heightΔhΔh=h∗−h0 being 0 mm, 10 mm, 20 mm, 30 mm, and 40 mm, respectively. Based on these, the initial porosity φ0 can be found according to (2).The heights of the naturally placed broken rock samples mixed according to different Talbol power exponents in the cylinder of the permeameter are different. The greater the Talbol power exponent is, the higher the height of the naturally placed sample is. The initial porosities of samples after compaction are different due to different compaction degrees. The scheme takes the gradation (Talbol power exponent) and the compaction degree (initial porosity) as two major considerations. A total of 50 groups of tests are required for 10 gradations and 5 compaction degrees. The test duration of each group should not be shorter than 18,000 s. At the beginning of the test, the lost masses are collected at a high frequency according to the quantities of migration and wastage of fine particles. And at the middle and later stages, the dissipated masses are collected every 3600 s.


## 2.1. Experimental Method and Principle

The patented pumping station seepage method was employed in our experiment on the mass-loss-concurring water seepage in broken rocks. Figure1 shows the experiment apparatus and Figure 2 shows the structure of its core device: an open permeameter.Figure 1
Schematic of the experimental setup.Figure 2
Schematic and photo of the open permeameter.The permeameter filled with broken rock sample is first installed on the electronic material testing machine. Then the porosity of the rock sample,φ∗, at its natural placement state can be obtained according to the following formula:
(1)φ∗=V∗−VzV∗=1−Mπa2h∗ps,where V∗ is the volume of the broken rock sample at its natural placement state, VZ is the volume of the compacted rock sample with the same mass, M is the total mass of the broken rock sample, ρs is the mass density of the broken rock sample, a is the inner radius of the cylinder, and h∗ is the height of the rock sample at its natural placement state.Then, an axial load is applied on the broken rock sample using the electronic material testing machine so as to control the amount of compaction and measure the height of the loaded rock sample,h0. Putting h0 into the following formula can find the initial porosity before seepage:
(2)φ0=V0−VzV0=1−Mπa2h0ρs,where V0 is the volume of the rock sample after it is forced to move a certain displacement along the axial direction under loading and h0 is the corresponding height.After the rock sample is immersed in water and saturated for half an hour, the plunger pump is turned on to push water in the water tank 1 continuously flowing crossing the rock sample from the bottom of the permeameter along the pipelines. As the water continuously flows crossing the rock sample for a long time, original and secondary fine particles in the broken rock sample migrate under the entrainment of water seepage, overflowing into the tray through the overflow tank of the permeameter and passing the hose into the vibrosieve.In the experiment, the filtration method is adopted to collect lost fine particles in a timely and stepwisely manner. During the test, the lost fine particles are collected once at a regular timeΔtii=1,2,3,…,n, with a total collected times note. After drying, the particles with a mass of Δmii=1,2,3,…,n are weighed out from the lost fine particles. Obviously, the total lost mass in the rock sample during each time period is
(3)m~=∑Δmi,and the mass-loss rate within unit time and unit volume, qi, is
(4)qi=ΔmiV0Δti.In the seepage process, with the continuous migration and wastage of fine particles in the broken rock sample, the porosity of the rock sample will also continuously change over time. The current porosity of the broken rock sample at any time during the test is calculated from the lost mass as follows:(5)φi=1−M−∑Δmiπa2h0ρs,i=1,2,3,…,n,Obviously, the time-variation features of the mass-loss rateqi and the porosity φi of the broken rock sample can easily be found with the migration and loss of fine particles in the experiment. Through the pressure sensors and flow sensors installed inside the pipeline and their connected data acquisitors and computer, the pore water pressure and flow can be timely monitored and acquired.


## 2.2. Experimental Scheme

The collapse column body often contains up to 70% of clay minerals as mudstone [35]. Therefore, mudstone is taken as the test samples with measured average uniaxial compressive strength of 9.85 MPa.In order to explore the effects of both the particle sizes (gradation) and the compaction degree (initial porosity) on the time-varying characteristics of the mass-loss rate of the broken rocks in the long-term seepage, the following experimental scheme is designed.Before the experiment, all the broken mudstone particles are sieved into 8 groups with their particle sizes as 0–2.5 mm, 2.5–5 mm, 5–8 mm, 8–10 mm, 10–12 mm, 12–15 mm, 15–20 mm, and 20–25 mm and their gradations are calculated according to the following Talbol continuous grading formula:(6)Pd=MdM×100%=dDn×100%,where P is the percentage of the particles mass whose diameter is less than d, Md is the mass of broken rock sample with particle size no greater than d, M is the total mass of broken rock sample, d is the particle diameter of the current broken rock sample, D is the maximum particle size of broken rock sample, and n is Talbol power exponent.Taking into account all the human, financial, and material resources and other factors consumed in the experiment, the values of Talbol power exponent,n, for the test samples are assigned as 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0, respectively. Based on this, the results of other gradation samples can be obtained by numerical computation, so as to not only save the cost of the test but also obtain more calculated data than the tested results.In the test, the total mass of mudstone filled in each cylinder is 2000 g. Figure3 shows the mass distribution of experimental samples with 8 kinds of particle sizes in 10 gradations. Figure 4 shows the test samples mixed uniformly at n=0.4 and n=1.0 according to the Talbol continuous grading formula.Figure 3
Mass distributions of test samples with various particle sizes in 10 gradations.Figure 4
The test samples mixed according to Talbol continuous grading formula.
(a)
n=0.4
(b)
n=1.0In the experiment, the broken rock samples are compacted to 5 different heights, as shown as the compaction heightΔhΔh=h∗−h0 being 0 mm, 10 mm, 20 mm, 30 mm, and 40 mm, respectively. Based on these, the initial porosity φ0 can be found according to (2).The heights of the naturally placed broken rock samples mixed according to different Talbol power exponents in the cylinder of the permeameter are different. The greater the Talbol power exponent is, the higher the height of the naturally placed sample is. The initial porosities of samples after compaction are different due to different compaction degrees. The scheme takes the gradation (Talbol power exponent) and the compaction degree (initial porosity) as two major considerations. A total of 50 groups of tests are required for 10 gradations and 5 compaction degrees. The test duration of each group should not be shorter than 18,000 s. At the beginning of the test, the lost masses are collected at a high frequency according to the quantities of migration and wastage of fine particles. And at the middle and later stages, the dissipated masses are collected every 3600 s.


## 3. Experimental Results and Analysis

### 3.1. Time-Varying Characteristics of Mass-Loss Rate

Figure5 shows the time-varying characteristics of the mass-loss rate in broken mudstone. From the figures, it is clear that there are two kinds of mass-loss rate curves with time. The first ones correspond to the test samples with gradations of n=0.1–0.6, and they have a rapid change stage and a slow change stage. The second ones correspond to the test samples with gradations of n=0.7–1.0, and they show three different stages, namely, the initial gradual change stage, the rapid change stage, and the slow change stage.Figure 5
Time-varying characteristics of the mass-loss rate in broken mudstone samples.
(a)
n=0.1
(b)
n=0.2
(c)
n=0.3
(d)
n=0.4
(e)
n=0.5
(f)
n=0.6
(g)
n=0.7
(h)
n=0.8
(i)
n=0.9
(j)
n=1.0The test samples with gradations ofn=0.1–0.6 at the early seepage stage contain a huge amount of original and crushing-produced fine particles. Under the action of pore water pressure, the mass-loss rate rapidly reduces because there are more migratable particles within a short time, as shown in Figure 6(b). With the seepage time increasing, these original and crushing-produced migratable particles in the samples become fewer and fewer and most migratable fine particles are from the secondary fine particles produced by long-term water-rock interaction. This stage lasts for a longer time and has fewer losable particles. Thus, the mass-loss rate at the stage shows a slowly varying characteristic, as shown in Figure 6(c).Figure 6
Schematic of migration and dissipation of fine particles in the test samples with gradations ofn=0.1–0.6.
(a)
Before seepage
(b)
Rapid change stage
(c)
Slow change stageThe test samples with gradations ofn=0.7–1.0 contain more particles of large sizes at a natural placement. After compaction, the fine particles with small sizes effectively fill the voids between large particles, resulting in a relative dense porous structure and very scanty effective seepage channels inside the samples under pore water pressure. At the early seepage stage, the migration and dissipation of fine particles last a relatively long time, thus, presenting an initial gradual change stage, as shown in Figure 7(b). Under the continuous action of the pore water pressure, multiple regionally connected seepage channels are gradually formed inside the samples, resulting in more lost fine particles in a short period of time and a rapid change in mass-loss rate, as shown in Figure 7(c). As the seepage channels are fully connected and maturely developed, most migratable fine particles are from the secondary fine particles produced by the long-term interaction between seepage and rocks, resulting in a slow variation in the mass-loss rate, as shown in Figure 7(d).Figure 7
Schematic of migration and dissipation of fine particles in the test samples with gradations ofn=0.7–1.0.
(a)
Before seepage
(b)
Initial gradual change stage
(c)
Rapid change stage
(d)
Slow change stage


### 3.2. Time-Varying Characteristics of Porosity

In the seepage process, the porous structure of samples will change accordingly with the migration and dissipation of fine particles. Figure8 shows the variational characteristics of porosity with time. From the figures, it is obvious that changes of porosity become less dramatic with Talbol power exponent increasing. Corresponding to the time-varying curves of the mass-loss rate, the time-varying curves of the porosity also fall into two kinds with gradations of n=0.1–0.6 and gradations of n=0.7–1.0.Figure 8
Time-varying characteristics of the porosity in broken mudstone samples.
(a)
n=0.1
(b)
n=0.2
(c)
n=0.3
(d)
n=0.4
(e)
n=0.5
(f)
n=0.6
(g)
n=0.7
(h)
n=0.8
(i)
n=0.9
(j)
n=1.0Jointly from both Figures6(b) and 7(c), it can be seen that the time-dependent curves of the porosity at the early seepage stage are affected by the Talbol power exponent and the compaction degree. Acted by the pore water pressure, fine particles inside the samples with gradations of n=0.1–0.6 rapidly migrate and drastically run away while a large number of fine particles originally filled the samples with gradations of n=0.7–1.0 continuously migrate and run away, leaving a large number of voids or empty cavities in the matrix. The remaining particles in the samples under the action of seepage entrainment and self-gravity re-fill those voids, resulting in a reorganization of the matrix structure and a great change in the samples’ porous structure.Jointly from both Figures6(c) and 7(d), it is clear that after a large-scale particle reconstruction inside the samples, only the nonmigratable matrix particles remain in the samples and the porous structure inside the samples basically remains stable, forming the water seepage channels. The long-term dissolution, erosion, and abrasion of water seepage to the mudstone makes the remaining skeleton particles more smooth or rounder, and only a very small number of secondary fine particles migrate and dissipate with water seepage, resulting in the porosity slightly increasing. It is deducible that with longer water seepage scour, the porosity of the sample will gradually enlarge with the growth rate gradually decreasing and eventually tending to a certain steady value φstable.


### 3.3. Expression of Mass-Loss Rate

Analyzing the time-varying characteristics of mass-loss rate shown in Figure5 finds that the mass-loss rate is an exponential function of time. In the process of water seepage in the broken rock sample, because the migration and dissipation of fine particles in the matrix could lead to porosity changes with time, the mass-loss rate is also a function of the porosity. According to the experimental test, the relationship of the mass-loss rate to both time and porosity can be expressed as
(7)qt=C1φte−χ1t+C2φstable−φt1−e−χ2t,where φstable is the sable value of the porosity dependent upon the Talbol power exponent and the initial porosity and φt is the porosity of the sample at the moment t. The unit of coefficients C1 and C2 is the same as that of q, and the unit of power exponents χ1 and χ2 is s−1. C1 and C2 as well as χ1 and χ2 can be determined using the genetic algorithm.


## 3.1. Time-Varying Characteristics of Mass-Loss Rate

Figure5 shows the time-varying characteristics of the mass-loss rate in broken mudstone. From the figures, it is clear that there are two kinds of mass-loss rate curves with time. The first ones correspond to the test samples with gradations of n=0.1–0.6, and they have a rapid change stage and a slow change stage. The second ones correspond to the test samples with gradations of n=0.7–1.0, and they show three different stages, namely, the initial gradual change stage, the rapid change stage, and the slow change stage.Figure 5
Time-varying characteristics of the mass-loss rate in broken mudstone samples.
(a)
n=0.1
(b)
n=0.2
(c)
n=0.3
(d)
n=0.4
(e)
n=0.5
(f)
n=0.6
(g)
n=0.7
(h)
n=0.8
(i)
n=0.9
(j)
n=1.0The test samples with gradations ofn=0.1–0.6 at the early seepage stage contain a huge amount of original and crushing-produced fine particles. Under the action of pore water pressure, the mass-loss rate rapidly reduces because there are more migratable particles within a short time, as shown in Figure 6(b). With the seepage time increasing, these original and crushing-produced migratable particles in the samples become fewer and fewer and most migratable fine particles are from the secondary fine particles produced by long-term water-rock interaction. This stage lasts for a longer time and has fewer losable particles. Thus, the mass-loss rate at the stage shows a slowly varying characteristic, as shown in Figure 6(c).Figure 6
Schematic of migration and dissipation of fine particles in the test samples with gradations ofn=0.1–0.6.
(a)
Before seepage
(b)
Rapid change stage
(c)
Slow change stageThe test samples with gradations ofn=0.7–1.0 contain more particles of large sizes at a natural placement. After compaction, the fine particles with small sizes effectively fill the voids between large particles, resulting in a relative dense porous structure and very scanty effective seepage channels inside the samples under pore water pressure. At the early seepage stage, the migration and dissipation of fine particles last a relatively long time, thus, presenting an initial gradual change stage, as shown in Figure 7(b). Under the continuous action of the pore water pressure, multiple regionally connected seepage channels are gradually formed inside the samples, resulting in more lost fine particles in a short period of time and a rapid change in mass-loss rate, as shown in Figure 7(c). As the seepage channels are fully connected and maturely developed, most migratable fine particles are from the secondary fine particles produced by the long-term interaction between seepage and rocks, resulting in a slow variation in the mass-loss rate, as shown in Figure 7(d).Figure 7
Schematic of migration and dissipation of fine particles in the test samples with gradations ofn=0.7–1.0.
(a)
Before seepage
(b)
Initial gradual change stage
(c)
Rapid change stage
(d)
Slow change stage


## 3.2. Time-Varying Characteristics of Porosity

In the seepage process, the porous structure of samples will change accordingly with the migration and dissipation of fine particles. Figure8 shows the variational characteristics of porosity with time. From the figures, it is obvious that changes of porosity become less dramatic with Talbol power exponent increasing. Corresponding to the time-varying curves of the mass-loss rate, the time-varying curves of the porosity also fall into two kinds with gradations of n=0.1–0.6 and gradations of n=0.7–1.0.Figure 8
Time-varying characteristics of the porosity in broken mudstone samples.
(a)
n=0.1
(b)
n=0.2
(c)
n=0.3
(d)
n=0.4
(e)
n=0.5
(f)
n=0.6
(g)
n=0.7
(h)
n=0.8
(i)
n=0.9
(j)
n=1.0Jointly from both Figures6(b) and 7(c), it can be seen that the time-dependent curves of the porosity at the early seepage stage are affected by the Talbol power exponent and the compaction degree. Acted by the pore water pressure, fine particles inside the samples with gradations of n=0.1–0.6 rapidly migrate and drastically run away while a large number of fine particles originally filled the samples with gradations of n=0.7–1.0 continuously migrate and run away, leaving a large number of voids or empty cavities in the matrix. The remaining particles in the samples under the action of seepage entrainment and self-gravity re-fill those voids, resulting in a reorganization of the matrix structure and a great change in the samples’ porous structure.Jointly from both Figures6(c) and 7(d), it is clear that after a large-scale particle reconstruction inside the samples, only the nonmigratable matrix particles remain in the samples and the porous structure inside the samples basically remains stable, forming the water seepage channels. The long-term dissolution, erosion, and abrasion of water seepage to the mudstone makes the remaining skeleton particles more smooth or rounder, and only a very small number of secondary fine particles migrate and dissipate with water seepage, resulting in the porosity slightly increasing. It is deducible that with longer water seepage scour, the porosity of the sample will gradually enlarge with the growth rate gradually decreasing and eventually tending to a certain steady value φstable.


## 3.3. Expression of Mass-Loss Rate

Analyzing the time-varying characteristics of mass-loss rate shown in Figure5 finds that the mass-loss rate is an exponential function of time. In the process of water seepage in the broken rock sample, because the migration and dissipation of fine particles in the matrix could lead to porosity changes with time, the mass-loss rate is also a function of the porosity. According to the experimental test, the relationship of the mass-loss rate to both time and porosity can be expressed as
(7)qt=C1φte−χ1t+C2φstable−φt1−e−χ2t,where φstable is the sable value of the porosity dependent upon the Talbol power exponent and the initial porosity and φt is the porosity of the sample at the moment t. The unit of coefficients C1 and C2 is the same as that of q, and the unit of power exponents χ1 and χ2 is s−1. C1 and C2 as well as χ1 and χ2 can be determined using the genetic algorithm.


## 4. Genetic Algorithm

The determination of optimal coefficientsC1best, C2best, χ1best, and χ2best in (7), that is, the expression of the mass-loss rate obtained based on the experimental results is the key for accurately calculating the mass-loss rate of the broken rocks. Based on the experimental results, we established the time series of mass-loss rate and calculated the optimal coefficients, C1best and C2best, as well as χ1best and χ2best, of 50 groups of tests using FORTRAN language program. In the following, we illustrate the implementation steps of the genetic algorithm by taking the test with Tabol power exponent of 0.8 and the compaction degree of 40 mm as an example.


### 4.1. Setting the Time Series Data File of Lost Mass

Convert the data acquired from the test with the first column as time (s) and the second column as lost mass (g) and save the file as the TXT file for future use by the program.


### 4.2. Assigning the Relevant Parameters

The value of each parameter in this example is that the (1) initial heighth=1.26×10−1m and diameter d=1.00×10−1m; (2) acquisition period τ=1.0s and data length N=8; (3) initial population scale kgroup=80, crossover individual logarithm kcop=25, crossover probability pc=0.8, and mutation probability pm=0.2, and the crossover point and mutation point are randomly produced and the lengths of genes C1,C2, χ1, and χ2 are 6; (4) the maximum generation algebra T=31 and minimum fitness s=2000; and (5) range of coefficient valuation: C11≤C1≤C12, C21≤C2≤C22, χ11≤χ1≤χ12, and χ21≤χ2≤χ22, among which, set C11=1.23×103, C12=8.54×104, C21=3.31×102, C22=9.61×105, χ11=2.83×10‐5, χ12=5.47×10‐3, χ21=0.31×10‐7, and χ22=2.11×10‐5.


### 4.3. Generating the Time Series of Porosity and Lost Mass, Respectively

Read the time series data file of lost mass obtained in Section4.1 and obtain the time series of porosity after simple computation. Figure 9 shows the output time-varying curves.Figure 9
Time-varying curves of porosity and lost mass obtained using the genetic algorithm.


### 4.4. Generating the Initial Population

Table1 gives the randomly generated initial population with size of 80.Table 1
Initial population.
Number
Genotype
Phenotype
Fitness
C1
C2
χ1
χ2
1
000101011010001000010110
2.35E + 03
1.05E + 05
1.36E−04
3.56E−06
9.80E + 00
2
111111101011010001010100
2.65E + 04
1.73E + 05
2.77E−04
3.24E−06
9.38E + 00
3
111110001000000011000110
2.61E + 04
3.24E + 04
5.71E−05
1.04E−06
9.46E + 00
... ...∗
... ...
... ...
... ...
... ...
... ...
... ...
78
000011110010101111000101
1.51E + 03
2.01E + 05
7.49E−04
8.86E−07
6.67E + 00
79
111010101101010001001000
2.44E + 04
1.81E + 05
2.77E−04
1.36E−06
6.76E + 00
80
011011001100110110110101
1.15E + 04
4.84E + 04
8.59E−04
8.43E−06
6.71E + 00
∗Due to limited space, data with population order numbers of 4–77 are not listed.


### 4.5. Generating and Determining the Values ofC1best, C2best, χ1best, and χ2best

First, select 50 individuals with the matching right using the random walk method from the initial population listed in Table1 and randomly pair them together. Second, randomly generate the crossover point for each gene and conduct crossover operation according to the crossover probability of pc=0.8. Third, randomly generate the mutation point for 5 genes of each individual after crossover operation and carry out the mutation operation according to the mutation probability of pm=0.2. Forth, decode the individual after mutation operation and compute the corresponding fitness Fitn. Fifth, find the individual sequence number with the highest fitness kbest using the sorting method. If Fitnkbest<s and the generation count ig<T, continue to propagate until Fitnkbest>s or ig=T and stop propagation. Sixth, decode the individual with decoding number kbest to obtain the optimal genes C1best, C2best, χ1best, and χ2best.After two generations of propagation, the individual (gene) withFitnkbest>s=2000 is obtained. Tables 2 and 3 show the genotypes, phenotypes, and fitnesses of individuals in the 1st generation and 2nd generation, respectively. Obviously, the fitness of the 2nd individual of the 2nd generation population shown in Table 3 is 2243, greater than 2000. Thus, the 2nd individual (gene) is the optimal gene with its C1best=3.60×104, C2best=2.53×103, χ1best=5.29×10‐4, and χ2best=1.00×10‐7.Table 2
The 1st generation population.
Number
Genotype
Phenotype
Fitness
C1
C2
χ1
χ2
1
100000010000000000000010
1.36E + 04
6.45E + 04
1.00E− 05
4.14E−07
1.77E + 03
2
000000010000000000000010
2.65E + 02
6.45E + 04
1.00E−05
4.14E−07
1.77E + 03
3
000000010000000001000000
2.65E + 02
6.45E + 04
2.57E−05
1.00E−07
1.78E + 03
... ...∗
... ...
... ...
... ...
... ...
... ...
... ...
48
000000010000000100000010
2.65E + 02
6.45E + 04
7.29E−05
4.14E−07
1.77E + 03
49
000000010000000001010000
2.65E + 02
6.45E + 04
2.57E−05
2.61E−06
1.77E + 03
50
000000010000000000000010
2.65E + 02
6.45E + 04
1.00E−05
4.14E−07
1.77E + 03
∗Due to limited space, data with population order numbers of 4–47 are not listed.Table 3
The 2nd generation population.
Number
Genotype
Phenotype
Fitness
C1
C2
χ1
χ2
1
100000010000100000000010
1.36E + 04
6.45E + 04
5.13E−04
4.14E−07
1.77E + 03
2
000000010000100000000010
3.60E + 04
2.53E + 03
5.29E−04
1.00E−07
2.24E + 03
3
000000010000000000000010
2.65E + 02
6.45E + 04
1.00E−05
4.14E−07
1.77E + 03
... ...∗
... ...
... ...
... ...
... ...
... ...
... ...
48
000000010000000100000010
2.65E + 02
6.45E + 04
7.29E−05
4.14E−07
1.84E + 03
49
000000010000000000010010
2.65E + 02
6.45E + 04
1.00E−05
2.93E−06
1.77E + 03
50
001000000000100001000000
2.65E + 02
6.85E + 04
1.00E−05
4.14E−07
1.77E + 03
∗Due to limited space, data with population order numbers of 4–47 are not listed.


### 4.6. Computing the Mass-Loss Rate and Comparing the Experimental and Theoretical Values

Solving the algebra equation by putting the best genesC1best, C2best, χ1best, and χ2best into (7) finds the calculated time sequence qii=0,1,2,…,N of mass-loss rate and compares the sequence with its corresponding experimental value. Figure 10 shows the compared results.Figure 10
Comparison of the experimental and calculated values from the time-varying curves of mass-loss rate.UsingFitnkbest can inversely find the errors between its experimental and theoretical values as follows:
(8)Err=1Fitnkbest4=122434=3.95×10−14.From the above example, evidently, if the range of reference values is reasonable, using this algorithm to propagate only 1 to 2 generations could obtain a satisfying calculated result with its fitness of individuals far higher than 2000 and with theoretical mass-loss rate curves almost fully consistent with the experimental mass-loss rate curves.


## 4.1. Setting the Time Series Data File of Lost Mass

Convert the data acquired from the test with the first column as time (s) and the second column as lost mass (g) and save the file as the TXT file for future use by the program.


## 4.2. Assigning the Relevant Parameters

The value of each parameter in this example is that the (1) initial heighth=1.26×10−1m and diameter d=1.00×10−1m; (2) acquisition period τ=1.0s and data length N=8; (3) initial population scale kgroup=80, crossover individual logarithm kcop=25, crossover probability pc=0.8, and mutation probability pm=0.2, and the crossover point and mutation point are randomly produced and the lengths of genes C1,C2, χ1, and χ2 are 6; (4) the maximum generation algebra T=31 and minimum fitness s=2000; and (5) range of coefficient valuation: C11≤C1≤C12, C21≤C2≤C22, χ11≤χ1≤χ12, and χ21≤χ2≤χ22, among which, set C11=1.23×103, C12=8.54×104, C21=3.31×102, C22=9.61×105, χ11=2.83×10‐5, χ12=5.47×10‐3, χ21=0.31×10‐7, and χ22=2.11×10‐5.


## 4.3. Generating the Time Series of Porosity and Lost Mass, Respectively

Read the time series data file of lost mass obtained in Section4.1 and obtain the time series of porosity after simple computation. Figure 9 shows the output time-varying curves.Figure 9
Time-varying curves of porosity and lost mass obtained using the genetic algorithm.


## 4.4. Generating the Initial Population

Table1 gives the randomly generated initial population with size of 80.Table 1
Initial population.
Number
Genotype
Phenotype
Fitness
C1
C2
χ1
χ2
1
000101011010001000010110
2.35E + 03
1.05E + 05
1.36E−04
3.56E−06
9.80E + 00
2
111111101011010001010100
2.65E + 04
1.73E + 05
2.77E−04
3.24E−06
9.38E + 00
3
111110001000000011000110
2.61E + 04
3.24E + 04
5.71E−05
1.04E−06
9.46E + 00
... ...∗
... ...
... ...
... ...
... ...
... ...
... ...
78
000011110010101111000101
1.51E + 03
2.01E + 05
7.49E−04
8.86E−07
6.67E + 00
79
111010101101010001001000
2.44E + 04
1.81E + 05
2.77E−04
1.36E−06
6.76E + 00
80
011011001100110110110101
1.15E + 04
4.84E + 04
8.59E−04
8.43E−06
6.71E + 00
∗Due to limited space, data with population order numbers of 4–77 are not listed.


## 4.5. Generating and Determining the Values ofC1best, C2best, χ1best, and χ2best

First, select 50 individuals with the matching right using the random walk method from the initial population listed in Table1 and randomly pair them together. Second, randomly generate the crossover point for each gene and conduct crossover operation according to the crossover probability of pc=0.8. Third, randomly generate the mutation point for 5 genes of each individual after crossover operation and carry out the mutation operation according to the mutation probability of pm=0.2. Forth, decode the individual after mutation operation and compute the corresponding fitness Fitn. Fifth, find the individual sequence number with the highest fitness kbest using the sorting method. If Fitnkbest<s and the generation count ig<T, continue to propagate until Fitnkbest>s or ig=T and stop propagation. Sixth, decode the individual with decoding number kbest to obtain the optimal genes C1best, C2best, χ1best, and χ2best.After two generations of propagation, the individual (gene) withFitnkbest>s=2000 is obtained. Tables 2 and 3 show the genotypes, phenotypes, and fitnesses of individuals in the 1st generation and 2nd generation, respectively. Obviously, the fitness of the 2nd individual of the 2nd generation population shown in Table 3 is 2243, greater than 2000. Thus, the 2nd individual (gene) is the optimal gene with its C1best=3.60×104, C2best=2.53×103, χ1best=5.29×10‐4, and χ2best=1.00×10‐7.Table 2
The 1st generation population.
Number
Genotype
Phenotype
Fitness
C1
C2
χ1
χ2
1
100000010000000000000010
1.36E + 04
6.45E + 04
1.00E− 05
4.14E−07
1.77E + 03
2
000000010000000000000010
2.65E + 02
6.45E + 04
1.00E−05
4.14E−07
1.77E + 03
3
000000010000000001000000
2.65E + 02
6.45E + 04
2.57E−05
1.00E−07
1.78E + 03
... ...∗
... ...
... ...
... ...
... ...
... ...
... ...
48
000000010000000100000010
2.65E + 02
6.45E + 04
7.29E−05
4.14E−07
1.77E + 03
49
000000010000000001010000
2.65E + 02
6.45E + 04
2.57E−05
2.61E−06
1.77E + 03
50
000000010000000000000010
2.65E + 02
6.45E + 04
1.00E−05
4.14E−07
1.77E + 03
∗Due to limited space, data with population order numbers of 4–47 are not listed.Table 3
The 2nd generation population.
Number
Genotype
Phenotype
Fitness
C1
C2
χ1
χ2
1
100000010000100000000010
1.36E + 04
6.45E + 04
5.13E−04
4.14E−07
1.77E + 03
2
000000010000100000000010
3.60E + 04
2.53E + 03
5.29E−04
1.00E−07
2.24E + 03
3
000000010000000000000010
2.65E + 02
6.45E + 04
1.00E−05
4.14E−07
1.77E + 03
... ...∗
... ...
... ...
... ...
... ...
... ...
... ...
48
000000010000000100000010
2.65E + 02
6.45E + 04
7.29E−05
4.14E−07
1.84E + 03
49
000000010000000000010010
2.65E + 02
6.45E + 04
1.00E−05
2.93E−06
1.77E + 03
50
001000000000100001000000
2.65E + 02
6.85E + 04
1.00E−05
4.14E−07
1.77E + 03
∗Due to limited space, data with population order numbers of 4–47 are not listed.


## 4.6. Computing the Mass-Loss Rate and Comparing the Experimental and Theoretical Values

Solving the algebra equation by putting the best genesC1best, C2best, χ1best, and χ2best into (7) finds the calculated time sequence qii=0,1,2,…,N of mass-loss rate and compares the sequence with its corresponding experimental value. Figure 10 shows the compared results.Figure 10
Comparison of the experimental and calculated values from the time-varying curves of mass-loss rate.UsingFitnkbest can inversely find the errors between its experimental and theoretical values as follows:
(8)Err=1Fitnkbest4=122434=3.95×10−14.From the above example, evidently, if the range of reference values is reasonable, using this algorithm to propagate only 1 to 2 generations could obtain a satisfying calculated result with its fitness of individuals far higher than 2000 and with theoretical mass-loss rate curves almost fully consistent with the experimental mass-loss rate curves.


## 5. The Mass-Loss Rate of the Actual Rocks

The gradient of pressure on the sampleGpspec in the seepage experiment is much larger than that on the actual rocks Gprock. To generalize the experimental results in the lab to the actual rocks, it is necessary to convert the mass-loss rate in the experiment into the mass-loss rate of the actual rocks, that is, to revise (7). In other words, we introduce an acceleration factor f to study its effects on C1, C2, χ1, and χ2:
(9)f=GpspecGprock.According to the experimental results and the optimal calculated results using the genetic algorithm, clearly, the acceleration factor,f, has little effect on χ1 and χ2 and can be neglected in the numerical calculation. Thus, we focus on the impact of the acceleration factor on C1 and C2. Since the process of mass loss from actual rocks often lasts for several months, due to limited manpower, financial resources, and material resources, we can only establish the relationships of C1prot and C2prot of rocks to C1 and C2 of rock samples by analyzing the mechanism of mass loss. Extrapolating the experimental results can get the following relationship:
(10)C1prot=C1fν1,C2prot=C2fν2,where C1prot and C2prot are the coefficients in the expression of the mass-loss rate of the actual rocks, respectively; ν1 and ν2 are the power exponents of the acceleration factor corresponding to C1prot and C2prot, respectively, both of which depend on the Talbol power exponent, n, and the initial porosity, φ0. By referring to the results of seepage tests, the ranges of ν1 and ν2 can be determined as ν1=−1.2to−0.9 and ν2=−1.2to−0.9. Thus, we obtain the time-varying mass-loss rate of the actual rocks as follows:
(11)qt=C1fν1φte−χ1t+C2fν2φstable−φt1−e−χ2t.


## 6. Algorithm Design of Mass-Loss Rate and Case Analysis

### 6.1. Algorithm Design

The mass-loss rate is calculated using the Lagrange discrete element method. In this method, the continuous medium is discretized into Lagrange grid. In brief, the actual rocks are divided intoN units with equal or unequal length. These units are numbered i,i=1,2,…,N and have N+1 nodes with coordinates xj,j=1,2,3,…,N+1, as shown in Figure 11.Figure 11
Lagrange unit partition of the actual rocks.The calculation of mass-loss rate in this paper only involves the calculation of physical quantities of the nodes, which are denoted using the core sign, superscripts, and subscripts. Among them, the core sign is enclosed in parentheses “()”; the superscripts indicate the time and the subscripts indicate the positions of the nodes (ordinal number). For example, assuming that the time response of the mass-loss rate starts to be calculated att0=0 with the time step length of τ, the mass-loss rate at time ti=iτ,i=1,2,3,… and the node j,j=1,2,3,…,N+1 is completely expressed as qjti and the iterative format of the mass-loss rate of actual rocks at the node can be constructed according to (11), that is,
(12)qjti=C1fν1φjtie−χ1t+C2fν2φstable−φjti1−e−χ2t,j=1,2,…,N+1.


### 6.2. Definite Solution Conditions

A cylindrical broken rock of2r=10m in diameter and H=60m in height is taken to simulate the actual collapse column. Figure 11 shows its related coordinate system using the rock’s bottom as the origin and the vertical upward direction as the positive direction of the O-x-axis.Due to limited test samples, the numerical calculation can be used to calculate the scheme of any gradation. In this case study, we used samples with gradations different from the tested samples. For example, the Talbol power exponentn=0.68, the porosity of the actual rocks φp0=0.318 at the reference pressure p0=0, the water’s mass density ρlp0=0.1×104kg/m3, and the actual rocks’ bottom pressures pbase=1MPa and the top ptop=0, respectively. The other related parameters including the actual rocks’ pore compressibility cφ, the liquid’s compressibility cl, and the dynamic viscosity μ are set according to the demand for calculating the mass-loss rate response.As shown in Figure11, the actual rocks are evenly divided into 30 units, N=30, and the coordinates of the nodes are
(13)xj=j−1H30,j=1,2,…,31.with the boundary conditions of
(14)pj=1=pbase,pj=31=0,∂V∂xj=1=0,∂V∂xj=31=0,and the initial conditions of
(15)pt=0=pbase+j−130ptop−pbase,Vt=0=0,φt=0=φp01+cφpbase+j−130ptop−pbase,(16)ρlt=0=ρlp01+clpbase+j−130ptop−pbase.


### 6.3. Calculation Steps and Result Analysis

The calculation steps are as follows:Step 1.
The best power exponent valuesχ1best and χ2best of the mass-loss rate in 50 samples and the matrix ML with 50 rows are constructed, as shown in (17), whose columns 1 to 6 are the Talbol power exponent n, initial porosity φ0, coefficient C1, coefficient C2, power exponent χ1, and power exponent χ2, respectively.
(17)ML=0.10.2798.98×1047.68×1021.12×10−32.18×10−50.10.2261.73×1051.11×1021.19×10−32.34×10−50.10.1651.75×1054.05×1021.06×10−33.26×10−5⋮⋮⋮⋮⋮⋮1.00.3023.76×1027.00×1035.21×10−44.06×10−71.00.2536.09×1034.00×1045.25×10−44.14×10−71.00.1972.48×1034.50×1045.22×10−45.60×10−7.Step 2.
The relationships ofC1, C2, χ1, and χ2 to the initial porosity φ0 are fitted by the exponential function under different Talbol power exponents.
(18)C1=A1eb1φ0,C2=A2eb2φ0,χ1=A3eb3φ0,χ2=A4eb4φ0,(19)MLF=0.16.53×1057.38×1033.53×10−37.77×10−5−6.71−12.83−5.03−4.770.22.05×1055.28×1031.16×10−21.79×10−4−10.2−5.29−3.77−7.030.34.85×1051.21×1043.70×10−33.00×10−6−4.13−6.28−5.17−3.45⋮⋮⋮⋮⋮⋮⋮⋮⋮0.81.36×1061.43×1042.25×10−32.61×10−7−31.2−17.61−15.87−8.950.97.19×1063.75×1032.51×10−31.97×10−4−36.71−5.24−14.51−16.721.04.13×1066.09×1055.17×10−37.35×10−6−31.5−12.56−9.3211.68.Step 3.
The univariate three-point unequidistant Lagrange interpolation method is used to calculate the model of any gradation except for the test samples. In this case, the Talbol power exponentn=0.68 is taken as an example. The values of coefficients A1,A2,A3,andA4 and the powers b1,b2,b3,andb4 under this gradation are calculated and listed in Table 4.Step 4.
Taking the parameters in Table4 as the known reference values, considering the effect of the acceleration factor on the coefficients C1 and C2, the mass-loss rate values at 31 nodes, whose coordinate is xj=j−1H/30,j=1,2,…31, could be calculated by (12), as shown in Figure 12.Table 4
Values ofA1,A2,A3,andA4 as well as b1,b2,b3,andb4 at n=0.68.
A1 (gs−1m−3)
A2 (gs−1m−3)
A3 (s−1)
A4 (s−1)
b1
b2
b3
b4
2.51 × 107
1.22 × 105
2.39 × 10−2
2.41 × 10−4
−24.2
−9.31
−11.8
−9.41Figure 12
Spatiotemporal distribution of mass-loss rate of the actual rocks.Analysis of mass-loss rate at node 1 in Figure12 finds that within 0≤t≤50s, the mass-loss rate increases from 5.97 kgs−1 m−3 to 6.96 kgs−1 m−3 with an increment of 16.6%. After that, it starts to decrease, from 6.96 kgs−1 m−3 to 5.75 kgs−1 m−3 during 50≤t≤600s with a decrement of 17.4%, from 5.75 kgs−1 m−3 to 1.31 kgs−1 m−3 during 600≤t≤4800s with a decrement of 77.2%, from 1.31 kgs−1 m−3 to 0.369 kgs−1 m−3 during 4800≤t≤8400s with a decrement of 71.8%, and from 0.369 kgs−1 m−3 to almost zero as t≥8400s.Figure13 shows the distributions of the mass-loss rate along the height of the rocks at the time t=0s, t=20s, t=50s, t=600s, t=4800s and t=8400s.Figure 13
Spatial distribution of mass-loss rate at different time points.
(a)
t=0s
(b)
t=20s
(c)
t=50s
(d)
t=600s
(e)
t=4800s
(f)
t=8400sFrom Figure13, it is obvious that at node 1 (x=0), only original fine particles migrate upwards and no other particles supply, while at other nodes, in addition to the upward migration of the original fine particles, many particles are recruited from other nodes to continuously replenish, thus, resulting in the maximum net mass-loss rate at node 1 at t=0s. With the seepage continuing, when the porosity at node 1 is close to φstable=0.386, migratable fine particles remain few, so the mass-loss rate at node 1 is smaller than that at most other nodes.Analysis by combining Figure12 with Figure 13 finds that in both 0≤t≤50s and 50s≤t<8400s, fine particles migrate and drain in the irregular and regionally connected seepage channels, showing the characteristics of irregular distributions of mass-loss rate at various nodes at different moments, as shown in Figures 13(a)–13(f), but the general trend is that with the increase of seepage time, the mass-loss rate changes significantly at various nodes. Thus, during this period, the porous structure in the rocks undergoes significant redistribution. After the period, that is, t≥8400s, with the seepage time continuing, the mass-loss rate decreases incessantly to almost zero at various nodes, showing very small decrement in amplitude and very small change in porosity. Thus, within this period, the porous structure in the rocks redistributes slowly; the porosity is close to φstable=0.386 at various nodes and the porous structure remains basically stable.


### 6.4. Effect of Talbol Power Exponent on Mass-Loss Rate

Figure14 shows the time-varying curves of mass-loss rate of the rocks at the bottom of the model (at node 1) with Talbol power exponent taken separately at 0.12, 0.32, 0.62, 0.72, and 0.92, which are obtained by using the computing program for mass-loss rate in the dynamic seepage system of broken rocks under the same definite solution conditions. Figure 14(b) is the amplified view of Figure 14(a) in the first 300 s.Figure 14
Time-varying curves of the mass-loss rate at different Talbol power exponents.
(a)
(b)From Figure14, obviously, when the Talbol power exponents are less than 0.6, the time-varying curves of mass-loss rate have both rapid and slow change stages. When the exponents are greater than 0.6, the curves have three stages, namely, the initial gradual change stage, rapid change stage, and slow change stage. At the early seepage stage, the mass-loss rate decreases with Talbol power exponent increasing. After a long period of seepage, the mass-loss rate is weakly correlated with Talbol power exponent.


## 6.1. Algorithm Design

The mass-loss rate is calculated using the Lagrange discrete element method. In this method, the continuous medium is discretized into Lagrange grid. In brief, the actual rocks are divided intoN units with equal or unequal length. These units are numbered i,i=1,2,…,N and have N+1 nodes with coordinates xj,j=1,2,3,…,N+1, as shown in Figure 11.Figure 11
Lagrange unit partition of the actual rocks.The calculation of mass-loss rate in this paper only involves the calculation of physical quantities of the nodes, which are denoted using the core sign, superscripts, and subscripts. Among them, the core sign is enclosed in parentheses “()”; the superscripts indicate the time and the subscripts indicate the positions of the nodes (ordinal number). For example, assuming that the time response of the mass-loss rate starts to be calculated att0=0 with the time step length of τ, the mass-loss rate at time ti=iτ,i=1,2,3,… and the node j,j=1,2,3,…,N+1 is completely expressed as qjti and the iterative format of the mass-loss rate of actual rocks at the node can be constructed according to (11), that is,
(12)qjti=C1fν1φjtie−χ1t+C2fν2φstable−φjti1−e−χ2t,j=1,2,…,N+1.


## 6.2. Definite Solution Conditions

A cylindrical broken rock of2r=10m in diameter and H=60m in height is taken to simulate the actual collapse column. Figure 11 shows its related coordinate system using the rock’s bottom as the origin and the vertical upward direction as the positive direction of the O-x-axis.Due to limited test samples, the numerical calculation can be used to calculate the scheme of any gradation. In this case study, we used samples with gradations different from the tested samples. For example, the Talbol power exponentn=0.68, the porosity of the actual rocks φp0=0.318 at the reference pressure p0=0, the water’s mass density ρlp0=0.1×104kg/m3, and the actual rocks’ bottom pressures pbase=1MPa and the top ptop=0, respectively. The other related parameters including the actual rocks’ pore compressibility cφ, the liquid’s compressibility cl, and the dynamic viscosity μ are set according to the demand for calculating the mass-loss rate response.As shown in Figure11, the actual rocks are evenly divided into 30 units, N=30, and the coordinates of the nodes are
(13)xj=j−1H30,j=1,2,…,31.with the boundary conditions of
(14)pj=1=pbase,pj=31=0,∂V∂xj=1=0,∂V∂xj=31=0,and the initial conditions of
(15)pt=0=pbase+j−130ptop−pbase,Vt=0=0,φt=0=φp01+cφpbase+j−130ptop−pbase,(16)ρlt=0=ρlp01+clpbase+j−130ptop−pbase.


## 6.3. Calculation Steps and Result Analysis

The calculation steps are as follows:Step 1.
The best power exponent valuesχ1best and χ2best of the mass-loss rate in 50 samples and the matrix ML with 50 rows are constructed, as shown in (17), whose columns 1 to 6 are the Talbol power exponent n, initial porosity φ0, coefficient C1, coefficient C2, power exponent χ1, and power exponent χ2, respectively.
(17)ML=0.10.2798.98×1047.68×1021.12×10−32.18×10−50.10.2261.73×1051.11×1021.19×10−32.34×10−50.10.1651.75×1054.05×1021.06×10−33.26×10−5⋮⋮⋮⋮⋮⋮1.00.3023.76×1027.00×1035.21×10−44.06×10−71.00.2536.09×1034.00×1045.25×10−44.14×10−71.00.1972.48×1034.50×1045.22×10−45.60×10−7.Step 2.
The relationships ofC1, C2, χ1, and χ2 to the initial porosity φ0 are fitted by the exponential function under different Talbol power exponents.
(18)C1=A1eb1φ0,C2=A2eb2φ0,χ1=A3eb3φ0,χ2=A4eb4φ0,(19)MLF=0.16.53×1057.38×1033.53×10−37.77×10−5−6.71−12.83−5.03−4.770.22.05×1055.28×1031.16×10−21.79×10−4−10.2−5.29−3.77−7.030.34.85×1051.21×1043.70×10−33.00×10−6−4.13−6.28−5.17−3.45⋮⋮⋮⋮⋮⋮⋮⋮⋮0.81.36×1061.43×1042.25×10−32.61×10−7−31.2−17.61−15.87−8.950.97.19×1063.75×1032.51×10−31.97×10−4−36.71−5.24−14.51−16.721.04.13×1066.09×1055.17×10−37.35×10−6−31.5−12.56−9.3211.68.Step 3.
The univariate three-point unequidistant Lagrange interpolation method is used to calculate the model of any gradation except for the test samples. In this case, the Talbol power exponentn=0.68 is taken as an example. The values of coefficients A1,A2,A3,andA4 and the powers b1,b2,b3,andb4 under this gradation are calculated and listed in Table 4.Step 4.
Taking the parameters in Table4 as the known reference values, considering the effect of the acceleration factor on the coefficients C1 and C2, the mass-loss rate values at 31 nodes, whose coordinate is xj=j−1H/30,j=1,2,…31, could be calculated by (12), as shown in Figure 12.Table 4
Values ofA1,A2,A3,andA4 as well as b1,b2,b3,andb4 at n=0.68.
A1 (gs−1m−3)
A2 (gs−1m−3)
A3 (s−1)
A4 (s−1)
b1
b2
b3
b4
2.51 × 107
1.22 × 105
2.39 × 10−2
2.41 × 10−4
−24.2
−9.31
−11.8
−9.41Figure 12
Spatiotemporal distribution of mass-loss rate of the actual rocks.Analysis of mass-loss rate at node 1 in Figure12 finds that within 0≤t≤50s, the mass-loss rate increases from 5.97 kgs−1 m−3 to 6.96 kgs−1 m−3 with an increment of 16.6%. After that, it starts to decrease, from 6.96 kgs−1 m−3 to 5.75 kgs−1 m−3 during 50≤t≤600s with a decrement of 17.4%, from 5.75 kgs−1 m−3 to 1.31 kgs−1 m−3 during 600≤t≤4800s with a decrement of 77.2%, from 1.31 kgs−1 m−3 to 0.369 kgs−1 m−3 during 4800≤t≤8400s with a decrement of 71.8%, and from 0.369 kgs−1 m−3 to almost zero as t≥8400s.Figure13 shows the distributions of the mass-loss rate along the height of the rocks at the time t=0s, t=20s, t=50s, t=600s, t=4800s and t=8400s.Figure 13
Spatial distribution of mass-loss rate at different time points.
(a)
t=0s
(b)
t=20s
(c)
t=50s
(d)
t=600s
(e)
t=4800s
(f)
t=8400sFrom Figure13, it is obvious that at node 1 (x=0), only original fine particles migrate upwards and no other particles supply, while at other nodes, in addition to the upward migration of the original fine particles, many particles are recruited from other nodes to continuously replenish, thus, resulting in the maximum net mass-loss rate at node 1 at t=0s. With the seepage continuing, when the porosity at node 1 is close to φstable=0.386, migratable fine particles remain few, so the mass-loss rate at node 1 is smaller than that at most other nodes.Analysis by combining Figure12 with Figure 13 finds that in both 0≤t≤50s and 50s≤t<8400s, fine particles migrate and drain in the irregular and regionally connected seepage channels, showing the characteristics of irregular distributions of mass-loss rate at various nodes at different moments, as shown in Figures 13(a)–13(f), but the general trend is that with the increase of seepage time, the mass-loss rate changes significantly at various nodes. Thus, during this period, the porous structure in the rocks undergoes significant redistribution. After the period, that is, t≥8400s, with the seepage time continuing, the mass-loss rate decreases incessantly to almost zero at various nodes, showing very small decrement in amplitude and very small change in porosity. Thus, within this period, the porous structure in the rocks redistributes slowly; the porosity is close to φstable=0.386 at various nodes and the porous structure remains basically stable.


## 6.4. Effect of Talbol Power Exponent on Mass-Loss Rate

Figure14 shows the time-varying curves of mass-loss rate of the rocks at the bottom of the model (at node 1) with Talbol power exponent taken separately at 0.12, 0.32, 0.62, 0.72, and 0.92, which are obtained by using the computing program for mass-loss rate in the dynamic seepage system of broken rocks under the same definite solution conditions. Figure 14(b) is the amplified view of Figure 14(a) in the first 300 s.Figure 14
Time-varying curves of the mass-loss rate at different Talbol power exponents.
(a)
(b)From Figure14, obviously, when the Talbol power exponents are less than 0.6, the time-varying curves of mass-loss rate have both rapid and slow change stages. When the exponents are greater than 0.6, the curves have three stages, namely, the initial gradual change stage, rapid change stage, and slow change stage. At the early seepage stage, the mass-loss rate decreases with Talbol power exponent increasing. After a long period of seepage, the mass-loss rate is weakly correlated with Talbol power exponent.


## 7. Conclusions

The collapse column is widely developed in the North China, which has the characteristic of a complex geological structure. Once it is disturbed by mining, the migration and loss of fine particles during seepage in the broken rocks become the intrinsic cause for collapse column activation and water inrush. The variation characteristic of mass-loss rate is not only changing in time but also are affected by the porosity of the broken rocks. In order to study the water inrush mechanism resulting from the mass migration and loss, we researched the variation characteristics of mass-loss rate in dynamic seepage system of the broken rocks through experimental, theoretical, and numerical methods. From the fundamental research, we obtained the conclusions as follows.Basing on the tests, the mass-loss rate can be expressed as a time-dependent, exponential function with its coefficients related to the porosity. After an acceleration factor is introduced into the experimentally obtained time-varying mass-loss rate expression, the mass-loss rate of the actual collapse column model is numerically calculated. The spatiotemporal distribution characteristics of the mass-loss rate curves are significantly affected by gradation. When the Talbol power exponents are less than 0.6, the time-varying curves have two stages and when the exponents are greater than 0.6, the curves have three stages. At the early seepage stage, the mass-loss rate reduces greatly, the porosity enlarges continuously, and the porous structure in the rocks redistributes violently. After long-term seepage, the distribution behaviors of the mass-loss rate are basically unaffected by gradation, the mass-loss rate almost equals zero, the porosity gradually reaches a certain stable value, and the porous structure of the rocks remains stable.



---

*Source: 7137601-2018-07-15.xml*