The Influence of Hydrothermal Effect on Migration of Contaminants in Porous Media

This paper extends a one-dimensional model of adsorbed contaminants in saturated porous media to establish a three-dimensional case suitable for unsaturated porous media. The effects of hydrothermal influences on the distributions of temperature and moisture, evolution of the bulk strain and migration of nuclides are simulated and analysed, and the coupled thermo-hydro-mechanical governing equations determined by the author are presented. The results show that the saturation in the buffer layer gradually rises to 0.96 from the initial value of 0.5 after approximately 450 years, and increases slowly in the initial stage of the simulation. The nuclide concentration in the solidified body first increases and then decreases, reaching a maximum of 173.05 mol/m3 in approximately 110 years. Moreover, most of the radionuclides are released by buffer layer adsorption and the adsorbed concentration is almost 100 times the solute concentration. In addition, the calculation also shows that the nuclides penetrate the buffer and reach the surrounding rock approximately 400 years later.


Introduction
In recent years, environmental geotechnical engineering problems, such as thermal energy storage, geothermal resource development, nuclear waste disposal, heating pipeline design, and urban landfills, have received widespread attention at home and abroad. This topic involves changes in percolation, heat transfer, mass transfer, and porous media skeletons in porous media [1][2][3][4]. Chen [1] discussed the biochemical reaction, skeleton deformation, pore water migration, solute migration and pore gas migration mechanism of soil, and established a model considering the biochemical reaction-skeleton deformation-water vapor migration-solute migration coupling. The model control equations define and clarify the model parameters and test methods. Zhang [5] extended the thermo-hydro-stressmigration coupling model and a two-dimensional finite element program from to perform a threedimensional analysis. From the perspective of methodology research, a simple nuclear waste geological disposal model is taken as an example to carry out three-dimensional numerical simulation of the thermo-hydro-stress-migration coupling process, in which the near-field temperature, saturation and nuclide concentration are investigated to determine the distribution of and variation in the pore water pressure, displacement, normal stress, flow rate, etc. Feng et al. [6] provided basic research ideas for the study of thermo-hydro-stress-chemical coupling in the damage area of crystal rock excavation. The formation and evolution mechanisms of excavation damage zones are analysed, based on the experiments of typical underground laboratories (such as the Swedish Äspö hard rock laboratory). The elastic, elastoplastic, viscoelastic and plastic THMC analysis models and the development of efficient numerical analysis software systems were established. In addition, the thermo-hydro-stress-chemical coupling behaviour of the damaged area of crystal rock excavation was  [7] established a thermo-hydro-mechanical-chemical (THMC) coupled numerical analysis model to study the change of permeability of porous sedimentary rock composed of quartz, in which the effect of pressure on chemical processes was taken into account. Yin et al. [8] proposed an analytical model for the THMC coupling problem involved in CO 2 injection in the underground storage process based on the finite element method. Numerical experiments show that this model can successfully analyse the stress and pressure changes in the rock around the wellbore affected by thermal and chemical influences.
However, there are few studies on the quantitative analysis of the amount of adsorption that occurs during pollutant transport. Based on the established one-dimensional model of pollutant transport in saturated porous media, the governing equation suitable for three-dimensional unsaturated porous media is derived, and the expression of adsorbed concentration is given. In addition, the hydrothermal interaction effect on pollutant transport in unsaturated porous media can be effectively analysed by using the established governing equation of hydro-mechanical coupling in unsaturated porous media.

The Governing Equation
According to the one-dimensional model of pollutant transport considering adsorption in saturated porous media [9][10][11][12], the equation of pollutant transport considering adsorption in three-dimensional unsaturated porous media can be written as where C i is the solute concentration which represents the concentration of contaminants in the fluid (ML −3 ), C d is the adsorbed concentration which indicates the mass of the adsorbed solute in the solid phase per unit mass of porous medium(MM −1 ), u is the average fluid velocity in porous media(LT −1 ), D D is the mechanical diffusion coefficient(L 2 T −1 ); D e is the effective diffusion coefficient(L 2 T −1 ), S i is the source term (ML −3 ), ρ b is the bulk dry density (ML −3 ), and θ is the volumetric water content (LT -2 ). In Eq. (1), the mechanical diffusion coefficient tensor can be expressed as [13] ( ) where α L is the longitudinal dispersion (L), and α T is the lateral dispersion (L). The effective diffusion coefficient in unsaturated flow can be expressed as [14][15][16]  where D L is the diffusion coefficient of the pollutant (L 2 T −1 ), τ is the tortuosity of the pore, and n is the porosity of the porous medium.
Referring to the non-equilibrium adsorption model previously described [10][11][12], the nonlinear adsorption model [12] is used in this paper. Then, the solute concentration and the adsorbed concentration have the following relationship , 0 , 0 where k d is the adsorption equilibrium constant (L 3 M −1 ), β is the correction constant (L 3 M −1 ) of the equilibrium coefficient k NF , β≥0, and η is the proportional factor of the desorption equilibrium constant. When β = 0, Eq. (4) degenerates into a linear adsorption model.
The mass concentration and the molar concentration have the following conversion formulas.
where C i ' is the solute molar concentration (nL −3 ); C d ′ is the adsorption molar concentration which represents the number of moles of the solute adsorbed by the solid phase per unit volume of porous medium (nL −3 ), M N is the molar mass (Mn -1 ) . This paper assumes a molar mass of 1 g/mol. The impact of the process of pollutant transport on hydrothermal effects is not considered in this paper. Therefore, the hydrothermal coupling equation in unsaturated porous media [17] can be written as where u i is the displacement component (L), ε v is the Volumetric strain, p is the equivalent pore water pressure (ML -1 T -2 )(that is, p =s· p w +(1-s) · p g , p w is the pore water pressure (ML −1 T -2 ) and p g is the pore gas pressure (ML −1 T -2 )), s is the pore water saturation, T is the absolute temperature (θ) , ρ is the density (ML -3 ) ( the subscripts w, g, and s represent the pore water, pore gas, and porous medium matrix, respectively), θ is the volumetric water content, α B is the Biot coupling coefficient (L 2 T -2 ) and equals 1-K b /K m , β Tb is the thermal expansion coefficient of the whole porous medium (θ -1 ), β Mb is the wet expansion coefficient of the whole porous medium, λ and G are Lame constants (ML −1 T -2 ), K b is the bulk modulus of the porous media (ML −1 T -2 ) , K m is the bulk modulus of the solid matrix (ML −1 T -2 ), n is the porosity of the porous medium and g is the gravity acceleration (LT -2 ).
where θ sw is the saturated volume fraction of pore water, θ rw is the residual volume fraction of pore water, ξ α (L -1 ), ξ L , ξ n and ξ m are the coeffic ients of the VG model and ξ m =1-(1 /ξ n ).
where K is the permeability of the porous medium (L 2 ), K rw and K rg are the relative permeability of pore water and pore gas, μ is the fluid viscosity coefficient (ML -1 T -1 ), C p is the capillary action coefficient (M -1 LT 2 ) and D is equal to (0,0,1).
where e is the internal energy of a unit mass (L 2 T -2 ) and equals c v ·( T-T 0 ), h is the heat content of a unit mass (L 2 T -2 ) and equals c p · (T-T 0 ), T 0 is the initial absolute temperature, c v and c p are the constant volume and constant pressure specific heat capacities (L 2 T -2 θ -1 ), k t is the total thermal conductivity of the porous medium (MT -3 θ -1 ) and equals θ w ·k w +θ g ·k g +(1-n) · k s , v is the absolute speed (LT -1 ) ( the subscripts w, g, and s represent the pore water, pore gas, and porous medium matrix, respectively), β Tm is the thermal expansion coefficient of the porous medium solid matrix (θ -1 ), β Mm is the wet expansion coefficient of the porous medium solid matrix and Q is the total heat source strength (L -1 MT -3 ).

Calculation Model and its Parameter Selection
To simplify the calculation model and reduce the calculation amount of the model, and qualitatively analyse the migration process of pollutants, moisture and temperature under the condition of unbalanced layout conditions, considering the relevant literature [5,18], it is assumed that the multiple shielding system of high-level radioactive waste consists of glass solidified body, a buffer material and the surrounding rock. The size of the surrounding rock is 36 m x 20 m x 20 m. The glass solidified body and the buffer layer are cubic bodies whose edge lengths are equal to 1 m and 4 m, respectively. In this paper, three solidified solids were selected as the research object, so that the relative influence of the solidified solids on each other could be investigated comprehensively. The median distance of solidified solids to the radioactive waste was 8m. According to the symmetry, considering the X-Z plane as the symmetry plane, half of the model is taken as the calculation model. The plane view is shown in Figure 1.  The initial temperature of planar the surrounding rock, buffer layer and solidified body is 20 °C, the temperature above and below the surrounding rock (i.e., along planes Z=0 m and Z=20 m) is kept constant, equal to 20 °C, the solidified body releases heat with a constant power of 1000 W, and the thermal conductivity is 55.0 (W.m-1/K). The remaining boundary provides thermal insulation. The initial saturation of the surrounding rock and the buffer layer is 1 and 0.5, respectively. The boundary of the surrounding rock is a zero-flow boundary and cannot be deformed. The initial pollutant concentration in the surrounding rock, buffer layer and solidified body is 0 mol/m3 , and the pollutant flux around the boundary of the surrounding rock is 0. It is assumed that the solidified body releases contaminants at a strength of 3.75 × 10 -5 × 2 -t / 10000 a mol / (m 3 • s) (a represents one year). The model parameters are shown in Table 1, and the remaining parameters are detailed in the literature [5,12,17]. Using the COMSOL Multiphysics numerical analysis platform, the above model is  Table 1, in which the thermal expansion coefficient and wet expansion coefficient of solid matrix are the same as those of integral coefficient.  Figure 2 shows the temperature distribution at 1000 years. Comparing Fig. 2 and Fig. 6, it can be seen that the temperature distribution is similar to the nuclide concentration distribution. As shown in Fig.  2(a), the three-dimensional distribution is ellipsoidal, and the temperature of the solidified body is approximately 65 °C; the temperature decreases with distance from the heat source. The temperature equipotential surface at the top and bottom edges of the surrounding rock is relatively flat. Figure 2(b) shows that the temperature gradient in the X-direction is larger than the gradient in the Y-direction, the buffer layer temperature is approximately 40 °C, and the surrounding rock temperature is approximately 28 °C. The evolution of temperature over time at points A~E (see Figure 1) is shown in Figure 3. It can be seen from the figure that the temperature rises rapidly in the initial stage and reaches the maximum temperature in approximately 10 years. Then, the temperature gradually stabilised, and the maximum temperatures of A~E are 72.97 °C, 46.76 °C, 31.05 °C, 46.13 °C and 27.69 °C, respectively. The corresponding steady temperatures are 66.07 °C, 44.55 °C, 31.00 °C, 43.98 °C and 27.69 °C. In the evolution process, the temperature difference between points C and E is large, which is mainly reflected in the temperature rise. Thereafter, the temperature difference changed slightly and remained at approximately 3.3 °C. The main cause of the temperature difference is the interaction between the two solidified bodies. In addition, the temperature at points C and E rises to the highest temperature and remains unchanged, and there is no temperature reduction phase at points A, B, and D.  Figure 3. The variation in temperature with time at the inspection points. Figure 4 shows the evolution of saturation with time at the observation points B~E. It can be seen from the figure that the saturation in the buffer layer increases gradually from 0.5 to 0.96 after 450 years, while the saturation of the surrounding rock changes little. In the initial stage (t < 10 a), the saturation increases rapidly, and the pore water in the surrounding rock flows to the buffer layer. Moreover, because of the high pore water pressure gradient, the flow velocity is high. With the flow of pore water, the saturation in the buffer layer increases, and the change rate decreases gradually (10 a < T < 450 a). After 450 years, the saturation of the buffer layer is the same as that of the surrounding rock.  Figure 5 shows the evolution of body strain with time at points A~E. It can be seen that the volume change in the surrounding rock and the solidified body is small (at points A, C and E), while the buffer layer (points B and D) has a relatively large bulk strain. The buffer layer first undergoes compression deformation, reaching a peak of -1.21% at 60 years. This is mainly due to the large thermal stress caused by the increase in internal temperature; the internal matter expands outward, squeezing the external matter. Subsequent volume expansion occurred, and the initial state was recovered at approximately 400 years. Although expansion occurs again in the subsequent stage, the expansion rate is significantly reduced. The body strain tends to be stable at approximately 800 years, and the body strains at B and D are 0.12% and 0.21%, respectively. This is because the internal stress tends to equilibrate gradually as the heat and humidity propagate outward, while the external stress gradually increases under the effect of the heat and humidity.  Figure 6 shows the distribution of the nuclide concentration C i ' at 1000 years. At this time, the nuclide concentration in the solidified body is approximately 88 mol/cm 3 , and the nuclide concentration shows an ellipsoidal distribution centred on the solidified body, as shown in Fig. 6 (a). The nuclides are concentrated in the buffer layer near the solidified body, the concentration in the surrounding rock is relatively low, approximately 2.3 mol/cm 3 , and a cross-effect layer is formed between the two solidified bodies. As shown in Fig. 6(b), the nuclides penetrate farther in the Y-direction than in the Xdirection. Those in the Y-direction penetrate approximately 5 m. This shows that the additive effect of pollution sources in the Y-direction is more obvious.   The adsorbed concentration C d ' at each position in Fig. 8 was significantly larger than the solute concentration C i ', indicating that most of the nuclides released from the solidified body are adsorbed by the buffer layer and the surrounding rock, and it can be seen in Fig. 8 that the adsorbed concentration in the buffer layer is much higher than that in the surrounding rock. Referring to Fig. 7, from the data of points B and D, the adsorbed concentration in the buffer layer rises rapidly during the phase of increased nuclide release from the solidified body, and the adsorbed concentration at point B at this stage is higher than that at D. At the subsequent stage, the rate of increase in the adsorbed concentration is slowed, and the adsorbed concentration values of B and D are approximately the same. The data of point C and D in Figure 8 indicate that the nuclides penetrate the buffer layer and reach the surrounding rock at around 400 years. The rate of increase in the adsorbed concentration at C is faster than that of E, which is due to the superposition of nuclide migration between the two solidified bodies.

Conclusion
A one-dimensional model of the contaminant transport in an established saturated porous medium is expanded to a three-dimensional case suitable for an unsaturated porous medium, and the specific expression of the adsorbed concentration is given. The calculations show that the model can accurately analyse the effect of the hydrothermal effect on the migration process of pollutants in unsaturated porous media.
The numerical analysis shows that the temperature in the solidified body and the buffer layer first increases and then decreases with time and that the interaction between the solidified bodies is more obvious in the surrounding rock than in the buffer layer. As the heat and humidity propagate outward, the near-field stress tends to gradually equilibrate, while the far-field stress gradually increases under the effect of the heat and humidity. The nuclides mainly concentrate in the buffer layer, and most of are adsorbed; the adsorbed concentration is almost 100 times that of the solute concentration.