#### Abstract

In the traditional heap leaching of rare earth minerals, the top of the rare earth pile is covered with leaching liquid. This creates trouble for vegetation restoration carried out timely on the top of the pile. In order to solve this trouble, a novel pile structure into which leaching liquid is laterally injected is proposed for heap leaching of rare earth. In this study, a laboratory test is carried out to study the formation and distribution of preferential ﬂow paths for the rare earth pile under a horizontal liquid injection condition. Furthermore, numerical simulations based on a dual-permeability model are conducted to investigate the influence of the preferential ﬂow paths on the seepage characteristics in the rare earth pile. The results show that, under the horizontal liquid injection condition, the fine particles of the rare earth move away from the liquid injection end and also toward the lower part of the pile. The migration of the fine particles results in the formation and connection of macropore, thereby generating preferential ﬂow paths in the rare earth pile. The preferential ﬂow paths are mainly distributed in the lower part of the pile near the liquid injection end. This causes the fluid in the lower part of the pile to seep faster significantly than that in the upper part. Within the region where the preferential ﬂow paths develop, the seepage in the early stage of the horizontal liquid injection is dominated by preferential ﬂow. The preferential ﬂow is more significant at the locations farther away from the liquid injection end.

#### 1. Introduction

China is one of the richest countries in the world in rare earth resources. There are a wide range of rare earth minerals in China [1]. In particular, ion-adsorbed rare earth minerals that are rich in medium and heavy rare earth elements are widely explored in Jiangxi, Guangdong, Fujian, and other provinces. At present, there are mainly two leaching methods for extracting ion-adsorbed rare earth elements from mines: the in situ leaching method and the heap leaching method [2]. Compared with in situ leaching, heap leaching is more controllable and can reduce pollution of leaching liquid to surrounding environments [3]. In addition, heap leaching is more efficient than in situ leaching in the extraction rate for different grades of rare earth [3]. Because of this, the heap leaching method is also widely used, especially in some rare earth mines with complex geological conditions that do not meet the requirements of in situ leaching.

In the traditional heap leaching of rare earth minerals, the top of the rare earth pile is covered with leaching liquid. This leaching method delays the time for ecological restoration on the top of the pile and increases the difficulty of ecological restoration [4]. Therefore, there is an urgent need to find a pile structure and liquid injection scheme that is suitable for ecological restoration. In this regard, we propose a novel growable pile structure for heap leaching of rare earth, as shown in Figure 1. In this novel pile structure, rare earth is stacked horizontally step by step, and leaching liquid is laterally injected into the pile from the starting end. After the leaching of each level pile is completed, the stacking and leaching of the next level pile starts. Meanwhile, vegetation restoration can be carried out on the top of the current level pile because it is free from leaching liquid. When using this growable pile structure for heap leaching of rare earth, the time of ecological restoration on the top of the pile can be advanced, thereby reducing the adverse impacts of pile leaching on surrounding environments.

An important problem of mineral heap leaching is the formation of preferential ﬂow paths in the pile, which adversely affects the leaching efficiency and pile stability. McBride et al. [5] numerically investigated the preferential flow behaviour in unsaturated packed beds and heaps by utilizing a robust computational fluid dynamics framework that accounts for local preferential ﬂow paths in the heap leaching system. They found that the nonuniform and adverse flow due to the preferential ﬂow paths reduces the ore-liquid contact and hence leads to lower metal recoveries. Luo et al. [6] carried out an infiltration test in an ion-adsorbed rare earth mine with the goal to study the effect of the heterogeneity of leaching liquid on the leaching rate. The test results show that, with an increase in the preferential flow fraction, the leaching rate and leaching depth in the rare earth ore body decrease. In reality, the network of interconnecting void space and the distribution of preferential flow paths in a heap are often unknown. Ignoring the preferential flow will clearly lead to an overestimation in the leaching efficiency, especially in the early stage of leaching [5]. Kukemilks et al. [7] investigated the influence of preferential flow on hillslope stability by using a fully coupled, physically based hydrogeological model. Their investigations suggest that a local-scale increase of pore-water pressure due to the preferential flow reduces the slope stability. Shao et al. [8] also emphasized that the slope stability would be overestimated if preferential flow is not considered. Therefore, with regard to the novel growable pile structure for heap leaching of rare earth, it makes sense to understand the formation and distribution of preferential flow paths for the leaching efficiency and stability analysis.

Many studies have been carried out on the formation mechanism of preferential ﬂow paths and related seepage characteristics in soil. It is generally acknowledged that preferential ﬂow paths are channels with high permeability caused by change of pore structures [9–12]. Wu et al. [10] researched the pore structure evolution of ore granular media during leaching by using X-ray computed tomography. They deemed that the pore structure change is mainly due to the overall subsidence of ore particles under the hydraulic, gravitational, and chemical actions and the migration and accumulation of fine particles under the seepage action of leaching liquid. Zamani and Maini [11] explained that the particle transport and attachment can change the pore morphology and consequently the porosity of the porous medium. Sirivithayapakorn and Keller [12] studied the transport of colloids in saturated porous media, and it is found that the preferential ﬂow paths become more significant for larger colloids and greater pressure gradients. Many techniques such as breakthrough curves, dye tracing, and CT scanning are employed to quantify preferential flow in soil at different scales [13, 14]. These surveys show that preferential flow may accelerate or delay the movement of matter depending upon the position relationship between the matter and the preferential flow paths. In addition, some numerical models such as the dual-permeability model, the fracture network model, and the active region model are used to evaluate the seepage characteristics of preferential flow [15–18]. These numerical approaches seem to well capture the nonequilibrium transport characteristics of preferential flow. However, the current research mainly focuses on the conditions with in situ leaching under vertical rainfall infiltration or heap leaching under vertical leaching liquid injections. The formation of preferential ﬂow paths under a horizontal liquid injection attracts much less attention. The novel growable pile structure is characterized by a horizontal injection of leaching liquid at the starting end to facilitate ecological restoration on the top of the pile. Therefore, it is of significance to study the formation mechanism of preferential ﬂow paths and the seepage characteristics under a horizontal injection of liquid in the growable pile structure. It can provide a reference for the later blockage of the preferential ﬂow paths so as to achieve high-efficiency leaching and pile stability.

In this study, a small-scale test for a rare earth pile under a horizontal injection of liquid at the starting end is carried out in the laboratory. Through the laboratory test, the formation and distribution of preferential ﬂow paths are observed, and the related formation mechanism is revealed by analyzing the particle gradation. On the basis of this, numerical simulations of the laboratory test are conducted by using multiphysics simulation software COMSOL. It aims to investigate the influence of the preferential ﬂow paths on the seepage characteristics in the rare earth pile under the horizontal injection of liquid. Note that the chemical reactions between rare earth and leaching liquid are not considered in this study. The leaching liquid in the laboratory test and numerical simulations is replaced with water.

#### 2. The Growable Pile Structure for Heap Leaching of Rare Earth

In order to facilitate ecological restoration on the top of the pile during heap leaching of rare earth, we develop a growable pile structure for heap leaching, as shown in Figure 1. First, a heap leaching site needs to be built against a natural mountain, in which a vertical impermeable bedrock is required as the starting end of the heap. On the surface of the vertical bedrock, a pipe network is arranged for injecting leaching liquid into the pile. Then, the first-level pile is stacked next to the vertical bedrock. The compactness and permeability coefficient of the soil should meet the requirements of leaching efficiency and pile stability. After that, leaching liquid is horizontally injected into the pile to start the leaching of the first-level rare earth pile. When the leaching of the first-level pile is completed, the second-level pile is stacked and leached next to the first-level pile. In the meantime, vegetation restoration work is carried out on the top of the first-level pile. Following the same procedures, the pile continues to grow, and the rare earth leaching and vegetation restoration are progressing step by step simultaneously. During the stepwise stacking and leaching process, the liquid injection pressure of the pipe network needs to be increased to ensure the normal leaching of the piles on the later place. When a level of the pile undergoes a significant reduction in leaching efficiency, the entire stacking and leaching process ends at this position.

In order to obtain the geometric parameters of the growable pile structure, the pile stability under different slope heights, slope angles, and pile thicknesses is studied by using the limit equilibrium analysis method. The results show that, for each level pile, the pile thickness has little effect on the pile stability, while the slope height and slope angle significantly affect the pile stability. Through optimization calculations, it is recommended that the thickness of each level pile is set to 3 m, and the slope height is limited within 10 m. The slope angle of each level pile should be gradually reduced to improve the stability of the entire pile. Under a slope height of 8–10 m, the slope angle is recommended to be 60–55° in the first-level pile, 55–50° in the second-level pile, and 50–45° in the following level piles. Our calculations also show that, under the premise of maintaining pile stability, the pile can grow to 6–8 levels when increasing the liquid injection pressure at the starting end to ensure leaching efficiency.

For the growable pile structure for heap leaching of rare earth, this study focuses on the formation of preferential ﬂow paths under a horizontal injection of liquid and its influence on the seepage characteristics. The stepwise stacking and leaching process as well as vegetation restoration are not considered in the present study. Therefore, in the following laboratory test and numerical simulations, a whole pile structure without sublevels is developed. Notwithstanding its limitations, it is still helpful for clarifying the formation mechanism of preferential ﬂow paths under a horizontal injection of liquid.

#### 3. Laboratory Test of Preferential Flow Paths

##### 3.1. Model of the Rare Earth Pile under a Horizontal Liquid Injection

A small-scale test is carried out in a transparent glass groove for investigating the formation of preferential flow paths in a rare earth pile under a horizontal liquid injection. The rare earth sample is taken from the Lingbei ion-adsorbed rare earth mine in Dingnan County, Jiangxi Province, China. Before stacking, the soil sample is dried, sieved, and mixed thoroughly. The specific gravity of the rare earth sample is 2.587. The soil particles with a diameter less than 0.075 mm account for 19.0%, within 0.075–0.25 mm account for 32.1%, within 0.25–2 mm account for 23.2%, and greater than 2 mm account for 25.7%. These soil particles are stacked evenly in the glass groove to form a rare earth pile, as shown in Figure 2. During the stacking process, the pile is compacted in layers with a thickness of 2 cm in each layer. The void ratio of the soil in the pile is 1.0. The small-scale rare earth pile measures 0.25 m high, 0.4 m long at the top, and 0.65 m long at the bottom. The width of the pile is 0.28 m, equaling that of the glass groove. A water tank is installed at the starting end for injecting liquid into the pile. Between the pile and the water tank, a porous acrylic plate is inserted for separation. During the test, the water head in the tank remains constant at 0.25 m, which is the same height as the pile.

##### 3.2. Distribution of Preferential Flow Paths

To facilitate analysis, horizontal and vertical grid lines and ten observation points are drawn on the glass groove within the range of the pile section, as shown in Figure 3. The grid spacing is 0.05 m. After the liquid injection begins, a preferential flow path (PFP) develops near No. 10 observation point in the lower part of the pile. Two hours later, this preferential flow path extends to a position 0.15 m away from the liquid injection end, as shown in Figure 3(a). It continues to grow to the vicinity of No. 9 observation point after 20 hours of the liquid injection (see Figure 3(b)). Since then, this preferential ﬂow path stops its growth. In the later leaching process, there are no new preferential ﬂow paths appearing. Under the horizontal liquid injection condition, the preferential flow paths mainly develop in the lower part of the pile near the liquid injection end.

**(a)**

**(b)**

**(c)**

It is seen from Figure 3 that the wetting front in the lower part of the pile moves ahead of the wetting front in the upper part. This is partly because the injection pressure of water increases linearly from the top to the bottom under the horizontal injection of liquid. When the preferential ﬂow path arises, the wetting front at the bottom is only 0.05 m ahead of that at the top. At *T* = 2 h, the wetting front at the bottom leads by 0.15 m. When the preferential flow path extends to No. 9 observation point at *T* = 20 h, the lead of the wetting front at the bottom has reached 0.20 m. The formation and growth of the preferential ﬂow path also cause the liquid at the bottom to flow faster. This rapid liquid flow due to the preponderance ﬂow paths reduces the time of contact between the leaching liquid and the rare earth particles. It leads to insufficient chemical reactions of leaching in local locations, thereby reducing the leaching efficiency of the pile. Therefore, it is of significance to study the formation mechanism of the preferential ﬂow paths and the seepage characteristics of the preferential flow in rare earth piles. Note that when the preferential flow path stops its growth after *T* = 20 h, the lead of the wetting front at the bottom is then reduced. This also demonstrates that the preferential flow path accelerates the liquid flow in the lower part of the pile.

##### 3.3. Formation Mechanism of Preferential Flow Paths

This section studies the formation mechanism of the preferential flow paths under the horizontal liquid injection through analyzing the particle gradation of the rare earth pile after leaching. After two days of the liquid injection, the saturation line in the pile is found to be stable, and at this time, the leaching area of the pile has reached the maximum. Then, the liquid injection is stopped, and the pile is sliced into layers from the top to the bottom. At those ten observation points, a quality of rare earth particles is taken at each point by using a ring knife. After drying these soil particles, sieve tests are carried out to obtain the particle gradation of the rare earth particles at each point after leaching, as summarized in Table 1. In this study, the soil particles with a diameter greater than 0.25 mm are classified as coarse particles, and the particles that are smaller than 0.25 mm in diameter are classified as fine particles. The percentage of the particle weight occupied by the fine particles at each observation point is shown in Figure 4.

At Nos. 1–6 observation points located in the region without significant preferential ﬂow paths, the weight percentages of the fine particles after leaching gradually increase with an increase in the distance from the liquid injection end. However, before leaching, the particle gradation of the rare earth particles at each point is substantially identical. This comparison shows that, during the leaching process under the horizontal liquid injection, the fine particles move horizontally away from the liquid injection end due to the action of seepage force. The percentages of the fine particles at Nos. 4–6 observation points in the middle layer are higher than those of the points at the same horizontal distance in the upper layer. This means that, in addition to the horizontal movement, the fine particles also migrate to the regions with lower location potential energy. Note that the weight percentages of the fine particles at Nos. 7 and 8 observation points in the lower layer are lower than that at No. 4 observation point in the middle layer. This is because Nos. 7 and 8 observation points are located near the foot of the slope, and the fine particles at these two points flow outside the pile. It leads to a reduction in the weight percentages of the fine particles at these two points. Before leaching, the weight percentage occupied by the fine particles smaller than 0.25 mm is 51.1％. After leaching, the weight percentages of the fine particles smaller than 0.25 mm at different points are all less than 51.1%. This also indicates that the fine rare earth particles migrate away and flow outside the pile during the leaching.

Nos. 9 and 10 observation points are located in the region with a significant preferential ﬂow path PFP. The percentages of the fine particles at these two points are 28.93% and 34.42%, which are much lower than those at the surrounding points. This also shows that the fine particles originally located in the preferential ﬂow path are partially carried away and deposited in the surrounding locations. When the gaps between the coarse particles in the surrounding locations are filled by the fine particles to a certain extent, a reduction in the connected pores prevents further growth of the preferential ﬂow path. Moreover, as the distance from the liquid injection end increases, the kinetic energy of water gradually decreases due to energy loss. The action of water on the soil particles is thus reduced accordingly. This is another reason for the stop of the preferential ﬂow paths. Therefore, there are no penetrating preferential ﬂow paths formed in the pile.

From the above particle gradation analysis, it is concluded that, under the horizontal liquid injection, the fine rare earth particles migrate in the direction away from the liquid injection end and also toward the lower part of the pile. The movement of the fine particles results in a lack of filling between the coarse particle skeletons at local locations. Then, the unfilled pores are connected to form preferential ﬂow paths. These fine rare earth particles migrate to other locations around, where they strengthen the filling of the gaps between coarse particles, thereby preventing the preferential ﬂow paths from further extending to the surroundings.

#### 4. Numerical Simulations of the Seepage Characteristics under Preferential Flow

##### 4.1. Dual-Permeability Model

In Section 4, the influence of the preferential ﬂow paths on the seepage characteristics in the rare earth pile is numerically studied by using the multiphysics simulation software COMSOL. In COMSOL, the Richards equation is utilized to handle the fluid seepage in porous media. However, a single Richards equation can only deal with the seepage problem in the presence of soil matrix flow. When both soil matrix flow and preferential flow exist, the dual-permeability model based on the Richards equation needs to be used. In the dual-permeability model, it is assumed that the medium can be considered as a superposition of two different porous systems at the same position, as shown in Figure 5. One is the matrix domain of soil aggregate with weaker permeability, and the other one is the preferential flow domain of large pores with stronger permeability [19]. The two domains are both considered as homogeneous media with independent hydraulic and solute transport properties. The coupling between the two domains is realized through a soil water transfer term, which allows a dynamic exchange of water between the soil matrix domain and the preferential flow domain.

The dual-permeability model of unsaturated flow can be written in the following equations:where subscript *f* denotes the preferential flow domain, subscript *m* denotes the soil matrix domain, is the volumetric fraction of the pore space occupied by the respective ﬂow domain, *θ* is the volumetric soil water content, *K* is the unsaturated hydraulic conductivity, *h* is the soil water pressure head, *S* is the degree of saturation, and is the soil water transfer term, which is defined as the volume of ﬂuid moving from the preferential flow domain to the soil matrix domain per unit volume of soil per unit time.

In this study, the second-order soil water transfer term modified by Köhne et al. [20] is used. It is calculated as follows:where *β* is a dimensionless geometry coefficient related to the shape of soil aggregates, *a* is the average soil aggregate radius, *K*_{a} is the unsaturated hydraulic conductivity at the interface between the preferential flow and the soil matrix, and *h*_{i} is the initial pressure head. It is assumed that the initial pressure head is equal for the preferential flow domain and the soil matrix domain.

The analytic formulas of Van Genuchten [21] are frequently used to define the hydraulic properties *S*, *θ*, and *K*. These equations arewhere *θ*_{s} is the saturated water content that is equal to the porosity of the porous medium, *θ*_{r} is the residual water content, *K*_{s} is the saturated hydraulic conductivity, and *α*, *m*, and *n* are constants that specify a particular medium type, in which *m* = 1 − 1/*n*. These parameters need to be determined by laboratory penetration tests.

The composite hydraulic properties of the soil are dependent on the respective domain properties through summation rules [22]:where subscript *c* denotes the composite domain of the preferential flow domain and the soil matrix domain, and *q* is the area-weighted fluid flux density.

##### 4.2. Numerical Model and Hydraulic Properties

A two-dimensional numerical model of the rare earth pile in the laboratory test is developed by using the multiphysics software COMSOL, as shown in Figure 6. The numerical model measures 0.25 m high, 0.4 m long at the top, and 0.65 m long at the bottom, remaining consistent with the experimental rare earth pile. The two-dimensional model is discretized into a mesh of triangular elements. The element size is 0.01 m. The model contains a total of 2370 discrete elements. The development of the preferential flow paths inside the pile cannot be easily obtained. According to the observation on the side of the pile in the laboratory test, the preferential ﬂow paths mainly develop in the lower part of the pile near the liquid injection end. Therefore, only the lower right part of the numerical model is set as the coupling of the preferential ﬂow domain and the soil matrix domain. The rest is set as the soil matrix domain. In the dual-permeability domain, the volumetric fractions and are considered. Because the dried soil is used in the test, the initial pore water pressure *p*_{i} = −1.0 × 10^{5} Pa is adopted in the numerical study.

**(a)**

**(b)**

According to equation (3), the constants *α*, *m*, and *n* can be obtained through fitting the relationship between the saturation degree *S* and the soil water pressure head *h* (or matric suction Ψ). Liu et al. [23] measured the matrix suction of unsaturated red soil indirectly by using the filter paper method based on the water-gas exchange relationship between soil and filter paper. By using the measured data, they developed the relationship between the soil water content and the matric suction. In the present study, the same filter paper method as used by Liu et al. [23] is employed to test the matrix suction of the unsaturated rare earth. The soil samples with different void ratios (i.e., void ratio *e* = 0.7, 0.8, 0.9, and 1.0) are prepared in the test. According to the designed soil water content, the required amount of water is added and fully mixed with the prepared soil samples. In order to ensure the water exchange balance between the filter paper and the soil, the test time is not less than ten days. The measured data pairs of saturation degree and matric suction for the soil samples with different void ratios are plotted in Figure 7. Fitting these data pairs by using the power function in equation (3) yields constants *α* and *n*. Constant *m* can be obtained via the relationship *m* = 1 − 1/*n*. The hydraulic properties of the rare earth sample at the void ratio *e* = 1.0 are selected for the soil matrix domain, as listed in Table 2. At the void ratio *e* = 1.0, the constants *α* and *n* are 0.353 kPa^{−1} and 1.229, respectively. The void ratio *e* = 1.0 is converted into the porosity of 0.50. Therefore, the saturated water content *θ*_{s} = 0.50 is considered. According to the soil-water characteristic curve in Figure 7, when the matric suction Ψ tends to 10^{4} kPa, the saturation degree *S* approaches 0.2. The saturation degree of 0.2 can be considered as the residual saturation. Consequently, the residual water content *θ*_{r} = 0.10 can be obtained by multiplying the saturated water content and the residual saturation. Although only the hydraulic parameters under the void ratio *e* = 1.0 are used in the numerical modeling, the measured data pairs of saturation degree and matric suction at the other void ratios are still presented in Figure 7. The purpose of this is to show that the measurement result under the void ratio *e* = 1.0 is reliable through the comparison of the overall variation laws between the different void ratios.

According to the particle gradation tests conducted in Section 3.1, the ion-adsorbed rare earth used in the laboratory tests is classified into coarse-grained silty sand. Therefore, variable head permeability tests are carried out to determine the saturated hydraulic conductivity *K*_{s}, as presented in Figure 8. The saturated hydraulic conductivities at different void ratios follow the power function relation *K*_{s} = 12.197*e*^{5.817}. This relationship consists of the experimental results of reference [24], indicating the reliability of the test results in this study.

From the test results in Figures 7 and 8, the hydraulic properties of the soil matrix domain can be determined. For the preferential flow domain, the determination of these parameters is mostly experience dependent [25]. The saturated hydraulic conductivity of the preferential flow domain is determined through the relationship with the saturated hydraulic conductivity of the soil matrix domain. According to Kozeny–Carman equation, the saturated hydraulic conductivity *K*_{s} is given by [26]where *ρ* is the density of water, *η* is the dynamic viscosity coefficient of water, *φ* is the porosity of soil, *τ* is the tortuosity of the medium, *s* is the specific surface area of the porous medium, is the gravitational acceleration, and *c* is the Kozeny coefficient that is approximately equal to 2.0.

The soil matrix domain and the preferential flow domain are located in the same spatial position, and the fluids in the two domains are also the same. Therefore, the parameters *ρ*, *η*, , and *c* in equation (7) are equal for the soil matrix domain and the preferential flow domain. It is assumed that the specific surface area *s* in the two domains remains the same. Then,

The tortuosity *τ* can be defined as the ratio of an average length of microscopic flow paths to the length of the system in the direction of the macroscopic flux [27]. For the preferential flow, the microscopic flow path approximates the direction of the macroscopic flux. Therefore, in this study, the tortuosity *τ*_{f} = 1.0 is considered for the preferential flow domain. The porosity of the preferential flow domain *φ*_{f} also approaches 1.0. In this case, the saturated hydraulic conductivity of the preferential flow domain is thus obtained by

According to the study of Koponen et al. [26], the tortuosity of the soil matrix domain *τ*_{m} is a function of porosity *φ*_{m}where *φ*_{c} is the percolation threshold and *b* and *d* are fitting parameters. The tortuosity is a monotonously decreasing function of the porosity at *φ*_{m} > *φ*_{c} and begins to diverge at the threshold *φ*_{m} = *φ*_{c}. Koponen et al. [26] fitted the simulated tortuosity *τ*_{m} as a function of the porosity *φ*_{m} by using equation (10) and obtained *φ*_{c} = 0.33, *b* = 0.65 and *d* = 0.19. Substitution of these parameters and *K*_{sm} = 1.20 × 10^{−4} cm/s into equations (9) and (10) yields the saturated hydraulic conductivity of the preferential flow domain *K*_{sf} = 1.86 × 10^{−3} cm/s. Unlike the matrix domain, the hydraulic parameters of the preferential flow domain cannot be directly measured. In this situation, referring to the work in [7, 8], the saturated water content *θ*_{s}, residual water content *θ*_{r}, and the constants *α* and *n* of the preferential flow domain are determined, as shown in Table 2. The residual water content of the preferential flow domain is much smaller than that of the soil matrix domain due to its stronger permeability.

##### 4.3. Validation of Numerical Modeling

Two tests are employed to validate the numerical modeling through comparing the simulation results with the test results. One is the soil column test conducted by Köhne and Mohanty [28], which is used to verify the accuracy of the dual-permeability model with the second-order water exchange term. The other one is the rare earth pile test carried out in this study, which is used to verify the rationality of the used hydraulic properties. Köhne and Mohanty [28] performed the water migration test in a soil column with a cylindrical macropore. During the soil column preparation, a metal tube of 2.4 cm diameter was vertically attached to the flow divider and coarse sand was poured into it. This formed a 2.4 cm diameter macropore filled with sand at a relatively low bulk density. Conceptually, the loosely filled macropore represents a preferential flow path. The column was then packed with dry-sieved and graded sandy loam soil at a constant bulk density, which represents a matrix domain. After the preferential flow and soil matrix domains were established, the metal tube that separates the macropore and soil matrix was pulled out from the soil column. In the column, five time domain reflectometry (TDR) three-prong probes were inserted into the sandy loam soil at different depths to measure water content in the matrix domain, as shown in Figure 9. Two sets of TDR coil probes were placed inside the macropore at two depths to measure water content within the preferential flow.

The COMSOL modeling of the water flow processes in the soil column is carried out in this study by using the dual-permeability model. The model size, hydraulic properties, initial conditions, and boundary conditions remain the same as the soil column test. Figure 10 presents the measured and simulated water contents at 35 cm depth in the matrix and preferential flow domains. In the soil matrix domain, the water content remains almost constant at 0.35. However, for the preferential flow domain, the measured results of the TDR coil probes show that the water content decreases from 0.35 to 0.10 within the first 10 minutes and then remains almost constant. The initial drop and subsequent constancy are accurately simulated by the dual-permeability model. The satisfactory agreement between the measured and simulated water contents indicates the dual-permeability model is feasible in simulating the preferential flow.

Figure 11(a) presents the observed wetting front in the rare earth pile after the appearance of the preferential flow paths in the lower part. The water content of the rare earth at this moment obtained by the numerical simulation is shown in Figure 11(b). In the numerical results, the wetting front at the bottom is ahead of that at the top by 0.25 m, which agrees well with the test result of 0.20 m. The shapes of the wetting fronts also reach a good agreement between the numerical and testing results. This indicates that the used numerical model and hydraulic properties are appropriate for simulating the water flow processes in the rare earth pile.

**(a)**

**(b)**

##### 4.4. Seepage Characteristics under Preferential Flow

Figure 12 shows the volumetric soil water content of the rare earth pile at different times under the horizontal liquid injection. The numerical results presented in this figure correspond to the dual-permeability model. As a comparison, the numerical results of the single-domain flow model are also given in Figure 13. In the single-domain flow model, the volumetric fraction occupied by the soil matrix domain, , is equal to 1.0. Accordingly, the volumetric fraction occupied by the preferential flow domain, , is zero. Therefore, a single Richards equation is used to deal with the seepage problem in the single-domain flow model. The hydraulic parameters of the single-domain flow model are the same as those of the soil matrix domain listed in Table 2. When preferential flow paths do not exist (see Figure 13), the wetting fronts at the top and bottom of the pile move forward at almost the same velocity. However, when preferential flow paths develop in the lower part of the pile, the wetting front in the lower part moves faster than that in the upper part, as shown in Figure 12. This is because the water flow in the macropore domain is faster than that in the soil matrix domain, particularly in the initial stage (see Figure 12 at *T* = 2 h). After a long period of water exchange between the preferential flow domain and the soil matrix domain, the wetting fronts in the lower part of the pile for these two domains tend to coincide, as shown in Figure 12 at *T* = 10 h and *T* = 20 h. From the comparisons between Figures 12 and 13, it is seen that due to the formation and growth of the preferential ﬂow paths, the water flow becomes faster in the rare earth pile, especially in the lower part. This probably causes insufficient chemical reactions between the leaching liquid and the rare earth in local locations, thereby reducing the leaching efficiency of the pile. Furthermore, the macropores as a result of the preferential ﬂow are likely to cause local soil collapse inside the pile, which deteriorates the stability of the rare earth pile.

**(a)**

**(b)**

Figure 14 shows the soil water transfer term between the preferential flow domain and the soil matrix domain at Nos. 9 and 10 observation points. A positive value of represents the water transfer from the preferential flow domain to the soil matrix domain. Conversely, a negative indicates opposite soil water transfer. It is seen that, in the initial stage, water is transmitted from the preferential flow domain to the matrix domain at both No. 9 and No. 10 observation points. However, the soil water transfer term is greater at No. 9 observation point that is farther from the liquid injection end. Then, an equilibrium of water exchange is reached at No. 9 observation point. This indicates that the preferential flow dominates the water flow at locations far away from the liquid injection end. At No. 10 observation point near the liquid injection end, the water exchange mode then quickly turns to water transfer from the soil matrix domain to the preferential flow domain. Furthermore, the absolute value of the negative water transfer term is much greater than the positive value. This means that the soil matrix flow is dominant in the rare earth pile near the liquid injection end. The laboratory test results in Figure 4 show that the weight percentage of the fine particles after leaching at No. 9 observation point is lower than that at No. 10 point. At No. 9 observation point, the water transfer term from the preferential flow domain to the matrix domain is greater. This comparison shows that, in the region where the preferential flow predominates, the weight fraction occupied by the fine particles is smaller.

The above numerical results reveal the influence of the preferential ﬂow paths on the seepage characteristics of the rare earth pile under a horizontal liquid injection. Note that, in the numerical modeling, the soil medium of the rare earth pile is considered homogeneous. However, with regard to actual rare earth mine piles, the soil structure is complex. Many uncertainties may be involved in the formation and distribution of the preferential ﬂow paths for the actual piles. These issues under a horizontal liquid injection condition need to be further studied in the future. Notwithstanding the limitations of the numerical modeling, it is still helpful for understanding the water flow processes in the rare earth pile under a horizontal liquid injection condition.

#### 5. Conclusions

In order to facilitate vegetation restoration on the top of the rare earth pile during heap leaching, a growable pile structure that is characterized by the horizontal liquid injection is proposed in this study. The laboratory test and the numerical modeling based on the dual-permeability model are carried out to study the formation mechanism of preferential ﬂow paths and its influence on the seepage characteristics in the rare earth pile. The results show that, under the horizontal liquid injection condition, the fine particles of the rare earth migrate in the direction away from the liquid injection end and also toward the lower part of the pile. The migration of the fine particles results in preferential ﬂow paths developing in the lower part of the rare earth pile nears the liquid injection end. As a result, the fluid seepage in the lower part of the pile is significantly faster than that in the upper part. This reduces the contact time between the leaching liquid and the rare earth at local locations, probably causing insufficient leaching of the rare earth pile. Within the region where the preferential ﬂow paths develop, the fluid flow processes in the early stage of the horizontal liquid injection are dominated by preferential ﬂow. Furthermore, the degree of the preferential flow is higher at the locations farther away from the liquid injection end.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the National Key Research and Development Program of China (2019YFC0605005 and 2019YFC0605001) and the National Natural Science Foundation of China (52179102 and 51969015). The authors wish to express their thanks to all the supports.