**1. Introduction**

Currently, lithium-ion batteries (LIB) are one of the most common devices for energy storage [1]. Many works are devoted to understanding the lithium intercalation/deintercalation mechanisms into electrodes and to optimization of LIB components for better electrochemical characteristics (see reviews [2,3]). The most widely accepted electrochemical model of LIB is the so-called pseudo-two-dimensional (P2D) model [4,5], where the intercalation and deintercalation in electrode particles, ion transport in electrolyte and separator are described in terms of the normal di ffusion equations. The generalization to the anomalous di ffusion case is proposed in [6]. However, in several electrode materials, the mutual solubility of the lithiated and delithiated phases is very low, and two phases coexist in a wide range of the state of charge (SOC). Such behavior is typical for charging and discharging in cathode based on lithium iron phosphate and in anode materials based on lithium titanate [7]. For this reason, lithiation and delithiation of the cathode particles can be associated with reversible phase transition in quasi-binary system, where the intercalant concentration can vary in the range from 0 to 1 during charging or discharging. The corresponding theoretical approach is based on the phase-field theory [8–16] that is commonly employed in studies on phase transitions [17–19].

The fact that the electrodes consist of particles of various sizes is often neglected. For example, the classic single particle model (SPM) and P2D model assume the electrodes consisting of spherical particles of identical size. Nonetheless, electrodes are typically made of porous materials with particles of di fferent sizes and shapes. As a result, the medium in LIB is a highly heterogeneous system. Neglecting the particle size distribution (PSD) largely underestimates the capacity in the case of elevated C-rates [20], and the SPM approaches poorly describe the battery performance at higher current densities. PSD is an important factor in the degradation of batteries [21,22], and it may be altered during battery operation due to cracking or agglomeration of particles [21]. In [23], the multiple-particle model indicates that the PSD broadening may lead to higher values of volumetric capacity and energy density. On the other hand, this broadening can cause amplification of electrode polarization [23].

There are several papers (see e.g., [20,23–26]) where the particle size dispersion is considered within the di ffusion approach. In most cases, the actual PSD is described by finite number of particle groups [20,25,26]. The multiple-particle approach meets computational di fficulties related to the need to solve di ffusion equations for each particle bin [26].

In the present paper, we extend the phase-field model based on the Cahn-Hilliard equation to simulate lithium intercalation dynamics in a cathode with particles of distributed size. Due to PSD, during charging/discharging, number of active electrode particles and their interfacial flux density depend on time. The boundary condition with time-dependent flux density for the Cahn-Hilliard equation is formulated from the universal laws established for speed of front propagation and evolution of single-phase concentration in individual particles. Calculation of the time-dependentflux density is performed by the developed general approach that can be employed for analysis of intercalation process in particle ensemble with arbitrary PSD.

#### **2. Phase-Field Intercalation in a Single Spherical Particle**

We start with a simplified phase-field model for a single submicron particle under galvanostatic condition. Expecting to establish some general rules of phase-field intercalation in a spherical particle, we analyze the dynamics for wide ranges of particle size and C-rate. The process of intercalation is simulated in terms of the Cahn-Hilliard (CH) equation [14,16,18,19]:

$$\frac{\partial \mathcal{C}}{\partial t} = M \nabla^2 \left[ \frac{\partial f}{\partial \mathcal{C}} - \kappa \nabla^2 \mathcal{C} \right] \tag{1}$$

where *c* ≡ *c*(*r*, *t*) is concentration field of intercalating atoms that defines the local state of charge of a cathode particle, *M* is the mobility, *f*(*c*) is the free energy of mixing in a quasi-binary system, κ is the gradient energy coe fficient. In this study, the free energy of mixing *f*(*c*) of quasi-binary solid solution is considered in the regular solution approximation [17–19]:

$$f(c) = \Omega c(1 - c) + k\_B T [c \ln c + (1 - c) \ln(1 - c)]\_+$$

Ω is the interaction parameter, *T* is the temperature, and *kB* is the Boltzmann constant. The influence of elastic strain is neglected. Assuming spherical symmetry of cathode particles, we use the Laplace operator of the form: ∇2 = ∂2 ∂*r*2 + 2 *r*∂ ∂*r*.

The insertion and extraction of intercalating atoms are simulated at the temperature of 300 K. The mobility of intercalated atom *M* in a cathode particle is related to the di ffusion coe fficient through the Einstein relation, *M* = *D*/(*kBT*). The di ffusion coe fficient of intercalating atom (e.g., lithium) for di fferent cathode materials is *D* = 10−13–10−15 m2/s at the considered temperature [4,8,13]. The interaction parameter Ω for di fferent cathode materials lies in the range of 0.059–0.193 eV/atom [8–16]. The gradient energy coe fficient can be related to the width of equilibrium profile [17,18] that is usually of several nanometers [8,10]. In this study, we employ the di ffusion coe fficient *D* = 10−14 m2/s, interaction parameter Ω = 0.115 eV/atom. The gradient energy coe fficient is usually written as κ = ηΩ*r*2 0 [8,14,17,18], where *r*0 is the intermolecular distance and η is the dimensionless coe fficient depending on form of interaction potential [18]. In this study, the value of gradient energy coe fficient is κ = 0.228 eVnm2.

The interaction parameter allows us to calculate phase diagram of the system (Figure 1). The equilibrium values of intercalating atom concentration are *cb*1 = 0.013 and *cb*2 = 0.987 at the temperature of 300 K. The unstable region lies in the interval between *cs*1 = 0.129 and *cs*2 = 0.871.

**Figure 1.** Phase diagram for regular binary solution with interaction parameter of Ω = 0.115 eV/atom.

The CH Equation (1) is solved under natural boundary conditions corresponding to the galvanostatic mode:

$$\left. \frac{\partial \mathbf{c}}{\partial r} \right|\_{0, \mathbb{R}\_0} = 0, \left. \frac{\partial \mu}{\partial r} \right|\_0 = 0, \left. \frac{\partial \mu}{\partial r} \right|\_{\mathbb{R}\_0} = -\frac{j\_r}{M}. \tag{2}$$

Here, μ = ∂ *f*/∂*c* − κ∇2*c* is the chemical potential defining the flux density of intercalating atoms. Initial distribution of the atoms is uniform and is characterized by equilibrium values of concentration, i.e., *c*0 = *cb*1 for insertion and *c*0 = *cb*2 for extraction process. The CH Equation (1) with boundary conditions (2) is solved numerically by using the explicit Euler time integration scheme [27].

Figure 2 demonstrates the evolution of concentration profile of intercalating atoms computed with the introduced values of the model parameter. The insertion and extraction processes are simulated for the cathode particle with radius of *R*0 = 0.1 μm under C-rate of 1C and 10C. Also, we simulated the insertion and extraction dynamics for cathode particle with radius *R*0 = 0.2 μm under C-rate equal to *C*/2.

Analyzing the results of insertion and extraction of intercalating atoms in cathode particles, we can identify some important features of these processes. At the very beginning of the processes, the concentration of intercalant decreases uniformly over the particle volume (see lines 1 and 2 in Figure 2a–f). At this stage, the particle corresponds to the single-phase pattern. The uniformity is explained by fast equilibration of atom distribution due to high diffusion coefficient and small size of the particle. The stage continues until the nearest value of metastability limit (*cs*1 or *cs*2) is achieved. At this stage, the average intercalant concentration *c* changes linearly with time (Figure 3). The approximate dependence can be easily obtained from the CH Equation (1):

$$
\overline{x} = c\_0 - \frac{3j\_r t}{R\_0} \,\text{\,\, \,}\tag{3}
$$

that agrees well with the results of direct solution of CH Equation (Figure 3). The relation (3) is valid at the time interval 0 ≤ *t* ≤ *ts*. From Equation (3) we can calculate the duration of the single-phase stage as *ts* = (*cb*2 − *cs*2)*R*0/3*jr* for extraction and *ts* = (*cs*1 − *cb*1)*R*0/3*jr* for insertion process.

**Figure 2.** Concentration profiles of intercalating atoms for different time points at *T* = 300 K. Corresponding radius of particle and *C*-rate are: (**a**) 0.1 μm, 1C insertion, (**b**) 0.1 μm, 1C extraction, (**c**) 0.1 μm, 10C insertion, (**d**) 0.1 μm, 10C extraction, (**e**) 0.2 μm, 0.5C insertion, (**f**) 0.2 μm, 0.5C extraction. The lines 3, 4, 9 and 10 correspond to passing from single-phase stage to double-phase stage and vice versa.

**Figure 3.** Kinetics of average concentration at the single-phase stage of insertion (1–3) and extraction (4–6) processes. Points are the solution of the CH equation, lines correspond to the linear dependence (3).

After that the concentration profile dramatically changes. Interfacial region depletes in very short time interval and the particle passes to double-phase pattern (see lines 3 and 4 in Figure 2a–f). At this stage, the intercalant redistributes over the particle volume and passes to the nearest equilibrium composition (see lines 4 and 5 in Figure 2a–f). Transition from single-phase to double-phase stage takes place for the short time interval estimated as ~0.2 s and ~0.5 s for particles with size of 0.1 μm (1C) and 0.2 μm, respectively. The longer time interval for larger particle is explained by the presence of additional concentration wave moving toward the particle center (see line 4 in Figure 2a–f) during equilibration of concentration profile. Overpotential of electrode depends on interfacial concentration of intercalating atoms [4,5,8]. Therefore, it is expected that passing from single-phase to double-phase regime can cause the abrupt change of overpotential. The effect could be used for determination of metastability limit.

The next stage is characterized by motion of the concentration wave from the particle interface to the center (lines 5–9 in Figure 2a–f). The position of concentration wave can be associated with coordinate *Rh*(*t*) corresponding to concentration *ch* = (*cb*1 + *cb*2)/2. The position of the wave front *Rh*(*t*) satisfies the linear dependence

$$\frac{R\_0^3 - R\_h^3(t)}{R\_0^2} = \frac{3|j\_r|t}{(c\_{b2} - c\_{b1})'} \tag{4}$$

obtained from the solution of the CH Equation (1). The Equation (4) is valid for *t* > *ts* (after the single-phase stage). The results of simulation by the CH Equation (Figure 2a–f) confirm Equation (4) (see Figure 4). The end of the double-phase stage corresponds to the time point *td* = *R*0(*cb*2 − *cb*1)/3*jr* for both insertion and extraction processes. The relation can be employed for calculation of the flux density *jr* corresponding to 1C-rate, if duration of charging or discharging process equals to *td* = 3600 s.

At the end of insertion or extraction process, the system passes over unstable region to another single-phase state that corresponds to fully charged or discharged state of a particle (lines 9 and 10 in Figure 2a–f). The process of insertion or extraction is stopped, when the corresponding equilibrium composition is achieved.

**Figure 4.** Kinetics of concentration wave propagation obtained by the CH Equation (points). Lines represent the linearized dependence (4). Size of particles and applied C-rate are indicated in Figure 2.

Contribution of the single-phase and double-phase stages can be estimated by ratio τ*s* = *ts*/*td*. This ratio can be easily found from Equations (3) and (4):

$$
\tau\_{s1} = \frac{c\_{s1} - c\_{b1}}{c\_{b2} - c\_{b1}}, \ \tau\_{s2} = \frac{c\_{b2} - c\_{s2}}{c\_{b2} - c\_{b1}}.
$$

for charging and discharging, respectively. The value of τ*s* is defined by the phase diagram and depends on temperature only. In case of symmetric phase diagram, this parameters are equal to each other τ*s*1 = τ*s*2. Calculation of the ratio in the range 250–400 K gives contribution of single-phase stage about 10–15%. Simulation of insertion and extraction for lower diffusion coefficient reveals deviation from concentration profiles at the single-phase stage only if the diffusion coefficient is less than ~10−17 m2/s.

Remarkably, variation of insertion and extraction rates in the range 1–10 C at the temperature of 300 K does not modify the two-stage mechanism and relations (3) and (4) remain valid (see Figures 2–4). Thus, Equations (3) and (4) can be used for determination of the composition profile at any time point of insertion or extraction process.

#### **3. Intercalation in Particles of Distributed Size**

Real cathodes consist of particles of various sizes and can be characterized by certain PSD. Under the assumption that the current is equally distributed over the active surface, smaller particles are charged or discharged for shorter time interval and withdraw earlier from the intercalation process. Larger particles have to balance the total current change in the galvanostatic mode. Therefore, flux density at the particle surface turns to depend on time *jr* = *jr*(*t*).

Let us modify the obtained relations with respect to size distribution of cathode particles. If the flux density at the particle interface depends on time after integration of CH equation we obtain for single-phase and double-phase regimes:

$$\begin{aligned} \overline{\boldsymbol{c}}(t) &= \mathbf{c}\_0 - \frac{3}{\mathbb{R}\_0} \int\_0^t \boldsymbol{j}\_r(t)dt, \ 0 \le t \le t\_s, \\\overline{\boldsymbol{R}}\_0^3 - \frac{\mathbf{c}\_k^3(t)}{\mathbf{c}\_0^2} &= \frac{3}{(\mathbf{c}\_{b2} - \mathbf{c}\_{b1})} \int\_0^t \left| \boldsymbol{j}\_r \right| dt, \ t\_s < t < t\_d. \end{aligned} \tag{5}$$

Whereas the total amount of intercalating atoms is conserved, the relation for total time of insertion or extraction process can be written in the form

$$\frac{1}{3}R\_0(c\_{b2} - c\_{b1}) = \int\_0^{t\_d} |j\_\prime(t)|dt.\tag{6}$$

Equation (6) can be used for definition of threshold value *R*min defining the minimal size of particle that can participate in intercalation process

$$\frac{d\mathcal{R}\_{\text{min}}}{dt} = \frac{\mathfrak{d}|j\_{\mathcal{I}}(t)|}{c\_{\mathfrak{k}2} - c\_{\mathfrak{k}1}}.\tag{7}$$

Then we assume that particles with *R* ≤ *R*min are excluded from the intercalation process (they are fully charged/discharged). The total current redistributes over the interface of residual particles with *R* > *R*min. Therefore, the galvanostatic mode can be determined by the following equation:

$$j\_r(t) \int\_{R\_{\min}}^{\infty} w(R) R^2 dR = j\_0 \{ R^2 \},\tag{8}$$

where *j*0 is the flux density at the beginning of intercalation process and *w*(*R*) is the size distribution function of cathode particles, *R*2 is the squared radius averaged over the whole ensemble of cathode particles.

Combining Equations (7) and (8) leads to the equation for the threshold value *R*min

$$\frac{d\mathcal{R}\_{\text{min}}}{dt} = \frac{3j\_0 \{\mathcal{R}^2\}}{\left(c\_{b2} - c\_{b1}\right) \int\_{\mathcal{R}\_{\text{min}}}^{\infty} w(R)R^2 d\mathcal{R}}.\tag{9}$$

Let us assume that ensemble of cathode particles is described by the gamma distribution

$$w(R) = \frac{R^{m-1}}{a^m \Gamma(m)} \exp\left(-\frac{R}{a}\right) \tag{10}$$

that is usually employed as PSD for different cathode materials [28,29]. Here *a* and *m* are the constant parameters and Γ(*m*) is the gamma function. The expectation value and dispersion of the distribution are *R* = *ma* and σ2 = *ma*2, respectively.

Integration of Equation (9) with respect to PSD (10) gives the implicit time-dependence of *R*min

$$\frac{a^3}{\Gamma(m)} \left[ \Gamma(m+3) + \frac{R\_{\text{min}}}{a} \Gamma \left( m + 2, \frac{R\_{\text{min}}}{a} \right) - \Gamma \left( m + 3, \frac{R\_{\text{min}}}{a} \right) \right] = \frac{3j\_0 \left\langle \mathbf{R}^2 \right\rangle t}{\left( \mathbf{c}\_{b2} - \mathbf{c}\_{b1} \right)}. \tag{11}$$

.

The threshold radius *R*min alters from zero to infinity, and the charging/discharging time is determined by

$$t\_{\max} = \frac{a^3 (c\_{b2} - c\_{b1}) \Gamma(m+3)}{3 j\_0 \langle \mathbb{R}^2 \rangle \Gamma(m)}$$

Figure 5 shows the time-dependence of threshold particle radius *R*min in dimensionless coordinates for different values of parameter *m*. The results in these coordinates do not depend on parameter *a* (see Equation (11)).

**Figure 5.** Dynamics of minimal radius *R*min of active particle in process of intercalation. Lines are the results of analytical solution (11) for different values of parameter *m*: (1) *m* = 2, (2) *m* = 4, (3) *m* = 8.

Analyzing Figure 5, we conclude that *R*min depends linearly on time at the beginning of insertion or extraction process that can be described by dependence

$$\frac{R\_{\text{min}}}{\Gamma(\text{R})} = \frac{\Gamma(m+3)}{m\Gamma(m+2)} \cdot \frac{t}{t\_{\text{max}}}.\tag{12}$$

The time interval corresponding to linear dependence is characterized by a constant flux density *jr*. After that the flux density *jr* of intercalating atoms grows very quickly and goes to infinity when time approaches to *t*max.

Taking Equation (7) into account, one can transform Equations (5) into

$$\overline{c}(t) = c\_0 - (c\_{b2} - c\_{b1})\frac{R\_{\text{min}}(t)}{R\_0}, \ 0 \le t < t\_{s\prime} \tag{13a}$$

$$\frac{R\_0^3 - R\_h^3(t)}{R\_0^2} = R\_{\min}(t), \quad t > t\_s. \tag{13b}$$

These equations describe the dynamics of intercalant distribution over the electrode particle at the single-phase and double-phase stage with respect to PSD.

Let us consider the impact of the PSD parameters on the extraction process from the particle ensemble characterized by PSD (10). Here, *R* = 0.1 μm and *m* = 4. We solve the CH Equation (1) for particle with size of *R*0 = 0.2 μm and time-dependent interfacial flux density *jr*(*t*) that defines the boundary condition (2).

The time-dependence of flux density *jr*(*t*) can be calculated by Equations (7) and (11). To represent *R*min(*t*) given by Equation (11) in explicit form, we choose the following approximation

$$\frac{R\_{\rm min}}{\langle R \rangle} = \sum\_{m=0}^{n} a\_{m} \left( \frac{t}{t\_{\rm max}} \right)^{2m+1} + b \, \text{arctanh} \left( \frac{t}{t\_{\rm max}} \right) \,. \tag{14}$$

where *n* is the order of approximation, *am* and *b* are fitting parameters. Using this approximation, we obtain the explicit dependence of the interfacial flux density for *n* = 0 in the form

$$j\_r(t) = j\_0 \left( 0.512826 + 0.487174 / (1 - (t/t\_{\text{max}})^2) \right).$$

The initial flux density *j*0 corresponds to C-rate of 1C for cathode particle with average size *R* .

The result of simulation of the extraction process from an individual particle of the ensemble is shown in Figure 6. The radius of this particle is chosen equal to *R*0 = 0.2 μm. The concentration profiles at different times are similar to Figure 2f. Analyzing Figure 6, one can conclude that solution of the CH Equation (1) with time-dependent interfacial flux density agrees well with the general Equation (13b) at the double-phase stage. The size distribution function of the cathode particles causes decrease of the total time of extraction that corresponds to time interval Δ*t* in Figure 6.

**Figure 6.** The time-dependence of concentration wave position *Rh* in dimensionless units for cathode particle with the size of *R*0 = 0.2 μm at the double-phase stage. The particle is the part of ensemble described by PSD given by Equation (10) with *R* = 0.1 μm and *m* = 4. Solid line 1 corresponds to solution of Equation (11). Dashed line 2 is the linear dependence corresponding to constancy of the interfacial flux density. Points 3 and 4 correspond to moments of start and finish of double-phase stage for the particle. Triangle points correspond to solution of the CH equation.

The time-dependent flux density is approximately constant at the single-phase stage, therefore the dependence of average concentration *c* at this stage is almost identical to that in Figure 2f. Generally, the time-dependent flux density can influence the average concentration *c* dynamics at the single-phase stage for large particles. The large particles are expected to deviate from the linear dependence that is observed under constant flux density (see Figure 3).


**1. Introduction**

Because of considerable migration depth, complex migration mechanism, sizeable spatial variability, strong concealment, and resistance to degradation of heavy metals (HMs), its high-concentration and high-risk pollutants in the soil are extremely difficult to repair and manage and are called "chemical time bombs" [1–3]. In recent years, the remediation of soil HM pollution has been an environmental topic worthy of attention. Studies have shown that HM pollutants may cause changes in soil structure, such as colloidal corrosion between soil particles, decrease in cementation strength, increase in clay content, resulting in increased soil compression, reduced liquid, and plastic limits, reduced shear strength (the reduction can reach about 35%), increased permeability, and the decreased bearing capacity [4–6]. In addition, the non-degradability of HMs leads to their accumulation in plants, animals, and humans along the biological chain [7–11]. Both pose a threat to human safety and health; therefore, it is imperative to control soil

**Citation:** Yang, Z.; Chang, J.; Li, X.; Zhang, K.; Wang, Y. The Effects of the Long-Term Freeze–Thaw Cycles on the Forms of Heavy Metals in Solidified/Stabilized Lead–Zinc–Cadmium Composite Heavy Metals Contaminated Soil. *Appl. Sci.* **2022**, *12*, 2934. https:// doi.org/10.3390/app12062934

Academic Editors: Bing Bai and Dibyendu Sarkar

Received: 6 February 2022 Accepted: 10 March 2022 Published: 13 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

HM pollution [12]. There are two kinds of solutions to manage HM pollution in the soil. The first one is to directly reduce the concentration of HMs in the soil, and the second one is to reduce the possibility of HM migration. That means fixing them in place and reducing their biotoxicity by special treatments [13,14]. In this context, a series of treatment methods for HM-contaminated soils have emerged, such as S/S technology, displacement technology, soil leaching, incineration, solvent extraction, glass curing, and physical separation, emerged. Among them, solidification/stabilization (S/S) technology is the most widely used because of its cost-effective and simple operation [14,15].

The S/S technology of soil usually provides the ability to handle HMs through two aspects. On the one hand, the S/S technology has a direct effect on HMs. Chemical precipitation, physical encapsulation, adsorption, and other means occur between HMs and hydrations products of binders; on the other hand, the S/S technology directly affects soil, and indirectly affects HMs by changing the soil structure. The reduced permeability leads to reduced migration of HMs, then the stability of HMs in the soil can be achieved [15,16]. However, while enjoying the restorative power of S/S technology, one must also consider its environmental impact.

Considering this problem, three inorganic binders, cement, lime, and fly ash, were used in the current study to improve S/S efficiency while reducing carbon emissions and easing the burden of industrial waste disposal [17–19]. However, it should be noted that the use of binders should involve reasonable control of the amount: too little S/S has poor effect, but if the amount is too much it will not only increase the cost, but also change the physical and chemical properties of the soil, and even cause secondary pollution.

Starting from the mechanism of the S/S technology, it can be realized that this technology cannot reduce the concentration of HMs in the soil, but restricts the migration of pollution, or limits the effectiveness of pollutants. Based on this, it is necessary to track the follow-up progress of HMs after the use of S/S technology. In the internal system of S/S HM-contaminated soil, the type of soil [20], the type and content of HMs [21], the type and content of binders [22,23], the curing conditions [23], and the moisture content [24] all have an impact on the long-term stability of HMs. In various complex environments, such as carbonation [25], acid rain leaching [25–27], high salt groundwater infiltration [28,29], wetting and drying cycles [30–32], and F–T cycles [33–36], the stable forms of the contaminants in the S/S-contaminated soil is very likely to change or even fail, which causes the activation of HM pollutants in the soil, causing secondary pollution and affecting its engineering mechanical properties. Therefore, the environmental durability of S/Scontaminated soil is extremely sensitive to environmental changes. The impact of F–T cycles on S/S HM-contaminated soil will continue to be focused on in-depth in this study.

Freeze–thaw cycles mainly occur in seasonal frozen regions. The moisture in the soil constantly transforms between liquid and solid with the periodic change of temperature: through changes in the volume of water in different fractions, as well as the free migration of water, the pore structure and density of the soil are changed, which is one of the direct factors affecting the migration effect of HMs. In China, the area of seasonal permafrost is about 5.14 million square kilometers, accounting for 53 percent of the total land area [37]. According to the second national soil erosion remote sensing survey, China's F–T erosion area has reached 13.4% of the country's total area [38]. Because of the vast F–T area, it is of great practical significance to explore the influence of the F–T cycle on the S/S of HM-contaminated soil. After 12 F–T cycles, the leaching concentrations of Pb/Zn-contaminated soil and composite Zn–Pb-contaminated soil increased greatly [34]. The elastic modulus of silty clay and modified clay can be significantly reduced after four F–T cycles [39]. After conducting seven F–T cycles on silty clay, Wang et al. concluded that the number of F–T cycles had a great influence on the mechanical properties of soil [40]. Studies have confirmed that under the long-term effect of simple F–T cycles (180 cycles, 180 days), the unconfined compressive strength loss rate of cemented contaminated soil can reach more than 30%, and the concentration of HM ions in HM leaching solution can increase by more than 20% [18,41,42].

This research has been improved on this basis, and "long-term" is the keyword of this. Studying the condition of HMs in the soil during the long F–T cycle period, observing the S/S effect of the HMs is a good way to ensure the long-term use of the S/S soil. Effective control of HMs in soil under adverse conditions is of great significance to prevent secondary pollution of HMs. A maximum of 90 F–T cycles (90 days) was set up under the same F–T mode to explore the stability of HMs under long-term F–T cycles in this research.

After making up for the problem of short-term F–T, the situation of HMs in the soil has also been improved, and the transformation from single HM pollution to compound HM pollution research has been realized. According to a large number of studies, one source of pollution often produces multiple types of HM pollution [41–43]. As, Cd, Cr, Cu, Mn, Ni, Pb, Zn, and other HMs exist in the farmland soil of the Qixia mining area in Nanjing, China [44]. There are many HMs in the weathered coal gangue, and the surrounding land is also harmed by HMs such as Zn, Cd, and Cr. Research on a single HM-contaminated soil can no longer meet the actual needs. Therefore, in this study, three pollutants, Pb, Zn, and Cd, were added to the soil at the same time.

The adsorption–desorption, precipitation reaction, complexation reaction, and oxidation–reduction reaction between HMs and soil allow HMs to exist in different forms in soil, and different forms of HMs exhibit different toxicity and environmental behavior [45]. When HMs are in a less bioavailable form or a more stable form, their risk in the environment is low. Therefore, the content of each form of HMs can be used as an observation indicator to assess the risk of HMs in soil. The improved Tessier method was used to extract seven fractions of HMs in the soil in this study.

Based on the above analysis, the improved Tessier method was adopted to detect the content of various forms of HMs (Pb, Zn, and Cd) in the soil that had undergone six F–T cycle conditions (0, 3, 7, 14, 30, and 90 cycles) up to 90 days, and the functional groups were analyzed by Fourier transform infrared spectroscopy analysis tests. Meanwhile, the study assumes that the HMs are uniformly distributed in the soil, by controlling the particles of the original soil less than 1 mm, the mixing time of the preparation of contaminated soil longer than 5 min, and the error of X-ray fluorescence analyses (XRF) of the same batch of three concentrations less than 3%.

In this study, the safety of the S/S technology is investigated by using a near-actual high concentration of complex hinged HM-contaminated soil, the morphological changes of HMs as a variable of interest, and the F–T cycle effects that may occur in most of the global land as the influencing factor. The conclusions and laws drawn from this study provide more rigorous theoretical support for predicting the safety of seasonally S/S HM-contaminated soil use and provide favorable conditions for post-remediation reuse of contaminated soil [18].