Thermo-Hydro-Mechanical Coupling Model of Unsaturated Soil Based on Modified VG Model and Numerical Analysis

A modified VG model considering both pore ratio and temperature effect is constructed. Based on the wet-thermal elasticity theory and mixture theory, coupled thermo-hydro-mechanical (THM) equations for moisture migration, heat transfer, and deformation in unsaturated soil are derived. The numerical implementation of the coupled theory was implemented by secondary development in the finite element platform. The validity of the established theoretical equation was verified by comparing it with the experimental data. Then the THM coupling response characteristics and internal mechanism of axisymmetric soil samples were numerically analyzed. Study shows that the water in unsaturated soil migrates from high-temperature zone to low-temperature zone under temperature load, and the permeability will change during the process of moisture absorption and dehydration. Particular attention should pay to the reasonable determination of the parameters in the modified VG model for the warming and cooling process of different measurement points during the numerical calculation. The higher the heat source temperature, the greater the saturation at the specific measurement point of the soil sample at thermal stability. The lower the heat source temperature, the shorter the time required for the saturation to stabilize at the specific measuring point. The volumetric strain of unsaturated soil results from a combination of wet expansion and thermal expansion, which is dominated by wet expansion near the heat source and mainly by thermal expansion far from the heat source. The change of the total displacement with time is dominated by the z-direction displacement, and its change trend along the radial direction is similar to that in the r-direction. The peak of the total displacement curve keeps moving backward with time.

A modified VG model considering both pore ratio and temperature effect is constructed. Based on the wet-thermal elasticity theory and mixture theory, coupled thermo-hydromechanical (THM) equations for moisture migration, heat transfer, and deformation in unsaturated soil are derived. The numerical implementation of the coupled theory was implemented by secondary development in the finite element platform. The validity of the established theoretical equation was verified by comparing it with the experimental data. Then the THM coupling response characteristics and internal mechanism of axisymmetric soil samples were numerically analyzed. Study shows that the water in unsaturated soil migrates from high-temperature zone to low-temperature zone under temperature load, and the permeability will change during the process of moisture absorption and dehydration. Particular attention should pay to the reasonable determination of the parameters in the modified VG model for the warming and cooling process of different measurement points during the numerical calculation. The higher the heat source temperature, the greater the saturation at the specific measurement point of the soil sample at thermal stability. The lower the heat source temperature, the shorter the time required for the saturation to stabilize at the specific measuring point. The volumetric strain of unsaturated soil results from a combination of wet expansion and thermal expansion, which is dominated by wet expansion near the heat source and mainly by thermal expansion far from the heat source. The change of the total displacement with time is dominated by the z-direction displacement, and its change trend along the radial direction is similar to that in the r-direction. The peak of the total displacement curve keeps moving backward with time.

INTRODUCTION
The THM multiphysics coupling process in unsaturated soil plays a significant role in engineering geology or environmental engineering . It has crucial applications in many areas, such as deep underground resource mining, geothermal resource development, nuclear waste disposal, and so on (Fu et al., 2020;Bai et al., 2021a;Bai et al., 2021b;Zhou et al., 2022). The theoretical model and numerical solution techniques for the THM coupling describing the characteristics of seepage, heat transfer, and mechanical behavior of unsaturated soils have become an important research content.
The research history of the THM coupling problem of unsaturated soil is relatively tortuous. The physical field develops from two fields to multiple fields, the influencing factors are considered from single to multiple, and the theoretical model is established from simple to complex. Booker (Booker and Savvidou, 1984) and Savvidou (Savvidou and Booker, 2010) theoretically solved the THM coupling problem using the mathematical transformation, and provided an analytical solution under the effect of point and cylindrical heat sources. Ghorbani et al. (2016) developed a theoretical model to completely describe the hydro-mechanical coupling phenomenon in non-saturation soil under static or dynamic loads, and verified the correctness of the established model through case analysis. Chen et al. (2009) established a theoretical model of the THM coupling problem in the three-phase system, and then the numerical analysis was performed based on the developed 3-D computer code THYME3D. Li et al. (2016) developed a mathematical model describing the THM coupling process by considering the three transport behaviors of convection, molecular diffusion, and mechanical escape of pollutants in unsaturated soil, and applied it to the simulation of an engineering clay barrier system. Akrouch et al. (2016) studied the heat exchange efficiency and hydrothermal coupling response characteristics of the heat source columns, and gave a calculation formula to quickly evaluate the heat exchange efficiency of heat source columns under different saturation levels. Recently, based on the theory of particle thermo-dynamics, Bai et al. (2021c) developed a new multi-field coupling theory considering the rearrangement of soil particles in the deformation process describing the microscopic changes of unsaturated soil particles under different loads. The soil-water characteristic curve (SWCC) is an effective mean of describing the hydraulic characteristics of unsaturated soil. It has been widely accepted in the theoretical derivation of the THM coupling problem by the numerous above scholars. However, the SWCC is difficult to obtain accurately due to multiple factors, resulting in incomplete numerical calculation and theoretical analysis.
Tarantinod (Tarantino, 2015) investigated the waterretention behavior of unsaturated soil across a wide scope of void ratios and then confirmed how void ratios affect the SWCC of unsaturated soil. Ghavam-Nasiri et al. (2019), Gallipoli et al. (2003), and Azizi et al. (2017) also analyzed the influence of void ratio on SWCC in different ways. The study reveals that the hydraulic properties can be more effectively described through establishing an SWCC that considers the pore ratio effect. In addition, temperature, as the main factor affecting the SWCC, has also been proved by many experiments. For example, Liu et al. (2002) conducted experimental research using an improved negative pressure soil tensiometer, indicating that temperature impacts the water potential of unsaturated soil. Cao et al. (2021) tested the SWCC of unsaturated soil under different heat loads and densities conditions using a filter paper-based method, and obtained the changing profile of SWCC. Romero et al. (2001), Jacinto et al. (2009), andLloret, (2004) carried out numerical experiments to study the effect of temperature on the moisture redistribution in seepage processes, and modified the classical SWCC considering the temperature effect to improve calculation accuracy. It is clear that the void ratio and temperature have a significant impact on the SWCC. However, it is rare to introduce the above two main factors simultaneously when modifying and optimizing the unsaturated SWCC.
Based on the classical VG model, this paper develops a modified VG model that takes the pore ratio and temperature effect into account. The THM coupling equation of unsaturated soil, which describes the moisture migration, heat transfer, and mechanical response in unsaturated soil, is derived based on the wet-thermal elasticity theory. The second developed COMSOL numerical platform was used to calculate the indoor tests to verify the constructed theoretical model. Then the coupled THM response characteristics of unsaturated soil was analyzed. van Genuchten, (1980) proposed a smooth and closed threeparameter mathematical model in 1980 to fit the SWCC. The expression is as follows:

MODIFIED VG MODEL
where α, n, and m are fitting parameters; S ew is the effective saturation; S r is the residual saturation; S w is the pore water saturation; S s is the maximum saturation; H c is matrix suction head; According to the work of Zhang et al. (2020), the SWCC fitting test of unsaturated soil was carried out based on Eq. 1. It can be observed from Figure 1 that the void ratio significantly affects the SWCC of unsaturated soil, and the increase in the void ratio causes the SWCC to shift. In addition, the SWCC curves for different pore ratios were plotted by fitting calculations using the experimental data (as shown in Figure 2). The effects of the three parameters α, n, and m on the shape variation of the VG model are given in Figure 3. According to Figure 3A, it can be seen the change of n causes SWCC to rotate around the turning point, and the change of m mainly affects the change of SWCC in the high suction region from Figure 3B. Comparing Figure 3C with Figure 2, it can be seen that there is a high similarity in the morphological changes of SWCC caused by the change of parameter α and pore ratio. In addition, according to the test of Tarantino, (2015), the relationship between void ratio and matrix suction approximately is a power function. Therefore, the correction considering the pore ratio effect on the classical VG model can start from the parameter α to introduce a correction term in the form of a power function. As opposed to establishing a connection between the parameters m and n (as shown in Eq. 2), the SWCC fitted by treating the parameters as independent variables agrees well with test data (Lu and Likos, 2004). The modified VG model considering the effect of void ratio is as follows: where e is void ratio; α m , n m, and m m are the fitting parameters of the modified VG model; b is the new proposed fitting parameter.
In numerical analysis, the influence of temperature on the VG model is generally reflected by fitting parameters α, n, and m FIGURE 2 | SWCC with different void ratio. The fitting calculations indicated that fitting parameters decrease linearly with the temperature increas (as illustrated in Figure 4). Therefore, it is assumed that the modified VG model parameters α m , n m , and m m are linearly related to temperature: where T is temperature; α 0 , n 0 , and m 0 are fitting parameters at the reference temperature T 0 ; λ α , λ n, and λ m are the fitting parameters of the modified VG model. Then the modified VG model considering the effects of temperature and void ratio can be obtained from Eqs 3, 4 as follows: where P c (P c = P g -P w ) is matrix suction; P w and P g are pore water and pore gas pressure, respectively; ρ w is pore water density; g is the acceleration of gravity.

Seepage Differential Equation
Unsaturated soil is a three-phase multi-component porous material made up of pore water, pore gas, and solid particles. According to the mixture theory (Malakpoor et al., 2007), The mass conservation law of each phase is as follows: where ϕ i , ρ i , and v i represent the volume fraction, density, and velocity of each component; The subscript i can be replaced by w, g, and s, representing pore water, pore gas, and solid particles, respectively. The volume fraction of each component is as follows: where n is porosity; S g is the saturation of pore gas, and S g = 1-S w . When pore water is incompressible and the gravity effect is considered, the velocity v of each phase in unsaturated soil has the following relationship (Rutqvist et al., 2001;Bai et al., 2020a): where K is permeability; K rw and K rg are the relative permeability of pore water and pore gas, respectively; μ w and μ g are the viscosity coefficients of pore water and pore gas, respectively; u w and u g are the average relative velocity of pore water and pore gas in the cross section of porous media, respectively. In addition, there are the following relationships (Qin, 2010): where ε v is volumetric strain. After the formula (11) is introduced into the mass conservation equation of solid particles, it can be calculated (Guo, 2018): Based on Eqs 8-12, the seepage differential equation of pore water for unsaturated soil can be obtained as follows: The seepage differential equation of pore gas is: where C p and C T are water-holding coefficients under capillary action and temperature action, respectively, which are two critical parameters of the seepage differential equation. They reflect the two driving effects in the seepage process and are determined by the SWCC (Guo and Bai, 2018). The following expressions can be obtained based on the modified VG model:

Heat Transfer Differential Equation
When ignoring the influence of soil particle displacement on heat transfer, the energy conservation equation can write as (Rutqvist et al., 2001;Bai et al., 2017;Bai et al., 2022): where K m is the bulk modulus of soil particles; e w , e g and e s are the internal energy per unit mass of pore water, pore gas and solid particles, respectively; β T and β M are the linear thermal expansion coefficient and linear moisture expansion coefficient of soil particles, respectively; h w , h g and h s are the heat content per unit mass of pore water, pore gas and solid particles, respectively; k t is the overall thermal conductivity; Q is the total heat source strength. The first item in Eq. 17 is the energy differential per unit time between entering and leaving; The second item is the energy caused by the fluid pressure change (Nowacki, 1970;Busch and Schade, 1976); The third item is the internal energy change rate in unsaturated soil; The fourth term is the change rate of reversible elastic energy due to strain of solid particles (Nowacki, 1970;Bear and Bachmat, 1990); The fifth item is the heat conduction term. The internal energy constitutive equation can be expressed as (Guo, 2018;Bai et al., 2020b): where e i is the internal energy per unit mass; h i is the heat content per unit mass; c vi and c pi are constant volume and constant pressure specific heat, respectively. The subscript i can be replaced by w, g, and s, representing pore water, pore gas, and solid particles, respectively. By substituting Eqs 9-12, 18, 19 into Eq. 17, the heat transfer differential equation can be obtained as follows: As a key thermodynamic parameter for numerical calculation, the overall thermal conductivity can be expressed as (Côté and Konrad, 2005): Where k s , k w , and k g are thermal conductivity of solid phase, liquid phase, and gas phase respectively; θ w is volumetric water content.

Principal Equation
According to the equivalent stress principle of unsaturated soil developed by Zienkiewicz et al. (1990), (Zienkiewicz et al., 1990), and Bear and Bachmat. (1990) and the effective stress increment equation developed by Rutqvist et al. (2001), the total stress in unsaturated soil is: where σ is the total stress; ε is total strain of solid skeleton; D and I are the tangential stiffness matrix and the unit matrix, respectively.‾P is weighted average pore pressure; α B is the Biot coefficient; ε P , ε T , and ε S are respectively the solid skeleton strains caused by pore pressure, temperature, and humidity, which are given by the following formula (Bai and Li, 2013;Bai et al., 2014;Fu et al., 2017;Guo, 2018): Then the stress-strain relationship in unsaturated soil is: Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 947335

NUMERICAL VERIFICATION AND COUPLING RESPONSE CHARACTERISTICS Physical Model and Parameters
To verify the reliability of the THM coupling theory constructed in this work, the numerical calculation was carried out based on the indoor test of unsaturated clay in the literature (Bai et al., 2020b), and the THM coupling response characteristics were analyzed. The physical model used in the calculation conforms to the axisymmetric (as given in Figure 5). The radius of the outer cylindrical surface of the sample is 0.5 m, and the height is 0.5 m; The radius of the heat source column is 0.02 m, and the height is 0.5 m. The symmetry axis of the heat source column coincides with the soil column. The four measuring points A, B, C, and D for temperature and moisture content were arranged sequentially outward from r = 0.078 m with a spacing of 0.1 m. The normal displacement of all boundaries of the hollow axisymmetric soil column (including the upper and lower base surfaces, inner and outer cylindrical surfaces) is zero. The soil under the action of the heat source is a closed system (no material exchange with the external environment). During the test, the top and bottom surfaces of the soil were adiabatic, and the outer cylindrical surface remained constant at 25°C. The initial temperature and initial volumetric water content of soil samples are 25°C and 19.6%, respectively. The temperature of the heat source rises to 80°C in a short time and remains constant for 5 days (a heating process), then falls to 25°C in a short time and lasts for 3 days (a cooling process). According to the relevant literature (Mualem, 1976;van Genuchten, 1980;Thomas and Vauclin, 1986;Guo, 2018), the parameters of the modified VG model are calculated as given in Table 1, 2 shows the parameters of the physical properties of the soil samples used in the experiments. Other state variables are calculated using a nonlinear relationship with temperature (Eckert & Rg, 1972;Vargaftik, 1975;Chen et al., 1999). Based on COMSOL numerical analysis platform, the corresponding program is compiled to model the THM coupling process in unsaturated soil. The seepage differential equation of gas adopts the Darcy seepage interface, the seepage differential equation of water and principal equation adopts porous elastic interface, and the energy equation adopts porous medium heat transfer interface. Due to the difference between the theoretical model and the built-in physical field equation in the COMSOL Multiphysics platform, the above physical field interfaces are secondarily developed by modifying the underlying equation to realize the numerical modeling of THM coupling theory. Figure 6 compares the numerical solution and the measured value of volumetric water content increment at different measuring points. The numerical results of measuring point A are lower than the measured data, and the numerical solutions of points B, C, and D are in close consistency with the experimental data. However, all simulation results are generally consistent with the variation trend of the volumetric water content incremental of the experimental data. Compared with the classical VG model, the numerical solution of the THM coupling equations using the modified VG model is in good agreement with the experimental data.

Model Validation and Analysis
Numerical simulation experiments show that the calculation results of volumetric water content increment in soil columns agree well with the experiment by changing the modified VG model parameters when simulating the corresponding cooling process of each measuring point. The specific approach is that the α 0 values of A, B, C, and D are 1.631, 1.623, 1.606, and 1.614, respectively. The volumetric water content of each measuring point grows continually during the heating phase of the soil sample (a hygroscopic process), and it gradually declines during the cooling process (a dehydration process). Because the SWCC of unsaturated soil during the moisture hygroscopic-dehydration process are different, i.e., the parameters of the modified VG model are different during the heating and cooling process of soil samples, and this difference must be fully considered in the numerical calculation process. In addition, it is found that the parameters of the modified VG model vary significantly when the distance between the measuring point and the heat source is different, and the closer to the heat source, the greater the difference in the model parameters. Figure 7 shows the comparison between the numerical results of temperature at different measuring points and the experimental values. It can be observed from Figure 7 that the calculation results of temperature at different points in all numerical calculations (both the modified VG model and the classical VG model) agree well with the experimental results. Overall, the calculated results of the temperature evolution of each measuring point during the heating process agree better with the experimental results than the cooling process. The deviation between the calculation results and the experimental results during the cooling process is significant, but the changing trend is consistent.
To quantify the difference between the predicted value of the theoretical model and the measured data, average relative error (AVRE) was used to evaluate the model's performance. The closer the AVRE value was to 1, the better the numerical simulation was. As can be seen from Table 3, compared with the classical VG model, the AVRE values of the numerical solution (including the volumetric water content and temperature during the heating and cooling processes) under the modified VG model are all closer to 1, indicating that the numerical simulation is better. In addition, the simulation values of temperature and volumetric water content during the heating and cooling of soil samples were compared, and the AVRE was calculated. The results show that the AVRE values of the heating process are closer to 1, indicating that the calculation accuracy in the process is higher than that in the cooling process. By comparing the AVRE values of temperature and moisture content under two kinds of SWCC models, it is found that the accuracy difference of temperature numerical solution is relatively small, indicating that the selection of the SWCC model has little effect on the heat transfer process. Figure 8 shows the saturation variation of measuring point B with time during the heating process of soil samples under different heat source temperatures. The temperature in unsaturated soil increases gradually under the continuous action of the heat source, which makes the water migrate from the high-temperature area to the low-temperature area, leading to a gradual increase in the saturation of point B. The increased rate of saturation of measuring point B is first fast and then slow. In Value 0.414 1 × 10 -22 1.5 2 × 10 -9 −1 × 10 -3 0.94 × 10 5 1.01 × 10 5 0.68  Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 947335 addition, the slope of the saturation curve at the initial stage of temperature increase is greater under higher temperature of heat sources. The higher the heat source temperature, the greater the slope of the saturation curve at the initial stage of the temperature increase. In the numerical calculation, the initial temperature of the soil sample is 25°C, and the set heat source temperatures are 35°C→50°C→65°C→80°C, respectively. Correspondingly, the time when the saturation of point B reaches a stable state is 0.2, 0.9, 1.3, and, and 2 days, respectively. It is known that the smaller the temperature difference between the initial temperature of the soil sample and the heat source, the faster the saturation of the soil sample reaches a stable state. Figure 9 shows the saturation contours of the symmetric plane at 0.01, 0.1, and 1 day during the heating process. It can be seen that in the initial stage (t = 0.01 days), the vertical distribution of soil saturation is slightly uneven due to the gravity effect in the range of r > 0.04 m. Still, there is no significant change relative to the initial saturation. When t = 0.1 day, the pore water moves outward from the vicinity of the heat source under the effect of the temperature gradient. The saturation of the soil sample in the left lower range (0.04 m < r < 0.35 m) increases gradually, showing a gradual decrease from the left lower to the right upper. When t = 1 day, the area where the saturation of the soil sample increases further expands, while the saturation near the heat source (r < 0.04 m) is still small. Figure 10A shows the distribution profile of saturation along the radial direction at z = 0.25 m at different times (namely, 0.1, 1, 1.5, and , 2 days) in the heating process of soil samples. It can be seen that the saturation curve moves upward with time, indicating that the saturation increases with time. When t = 1.5 days and t = 2 days, the saturation curve basically coincides, meaning that the soil saturation reaches a stable state when t = 2 days. In addition, the saturation curves generally show an accelerated rise followed by a slow decline, and the slope and the amplitude of change in the rising phase are larger. This phenomenon indicates that the saturation changes rapidly near the heat source, and reaches the maximum saturation at r = 0.02 m. The slope of the saturation curve becomes smaller as the radial distance increases in the descending section of the saturation curve. That is, the farther away from the heat source, the smaller the change of saturation along the radial direction, and the influence of the heat source on saturation becomes weaker with the increase of distance. Figure 10B shows the saturation variation of different measuring points with time during the heating process. It can be seen that the variation of saturation with time is similar to that of volumetric water content increment (as shown in Figure 6) with time. The saturation growth rate of measuring points A, B, C, and D showed the characteristics of first large and then small. The saturation of each measurement point to reach steady state was 0.544, 0.521, 0.506, and 0.492, and the increments of saturation were 0.074, 0.051, 0.036, and 0.022, respectively. It shows that the closer to the heat source, the greater the saturation increment during the heating process. Figure 11 shows the radial distribution of volumetric strain at z = 0.25 m at different times during the heating process of soil samples. The overall trend of the volumetric strain curve is highly consistent with the variation of saturation in the radial direction, and the maximum values of volumetric strain and saturation are in roughly the same position. The temperature gradient causes the pore water to migrate outwards and accumulates near the heat source to form a local high moisture content region. The wet swelling causes the volumetric strain of the soil in this region to increase gradually, and accordingly, volumetric strain reaches the maximum at the position of maximum saturation. In addition, the variation amplitude of volumetric strain in the decline stage is much larger than that of saturation, which is the combined effect of temperature and boundary constraints. The volumetric strain near the outer boundary of the soil column is negative, primarily passive compression deformation, and the volumetric strain in the middle of the soil is positive, primarily active expansion deformation. Figure 12 shows the distribution profile of temperature along the radial direction at different times during the heating process of soil samples. With the increase of time, the temperature gradually increases along the radial direction and reaches the thermal stability state approximately at t = 2 days. In addition, the temperature gradient near the heat source during  heat transfer decreases gradually with time. However, the temperature gradient away from the heat source increases gradually with time under the first type of thermal boundary. The comparison shows a certain similarity between the changes in temperature and volumetric strain along the radial direction in the region r > 0.3 m (That is, both are nearly straight at t = 0.1 day, while both curves at the remaining three times have a certain slope). The thermal expansion caused by the temperature difference in the region far from the heat source is the main factor driving the soil volumetric strain. Figures 13A-C show the radial distribution of total displacement, r-direction displacement, and z-direction displacement corresponding to four time points (0.1, 1, 1.5, and 2 days) in the heating process. It can be observed from Figure 13 that the total displacement increases first and then decreases along the radial direction, and decreases with time, which has the same variation profile with the displacement in the r-direction. The z-direction displacement is negative (the direction is opposite to the z coordinate direction), and the absolute value of the displacement increases along the radial direction and decreases with time. Because the absolute value of z-direction displacement is much larger than that in r-direction, the total displacement is mainly dominated by z-direction displacement. By comparing Figures 13A,B, it can be observed that the variation trend of the total displacement curve is similar to that of the r-direction displacement, and the peak points on the two displacements curves move outward with time.

CONCLUSION
A modified VG model considering both void ratio and temperature effect is proposed. The THM coupling theoretical model of unsaturated soil is given based on the wet-thermal elasticity theory and mixture theory. The numericalization of the theoretical model is realized by the secondary development of COMSOL Multiphysics. The validity of the theoretical model equation is verified using the indoor test data. Finally, the thermal-hydro-mechanical coupling response characteristics are analyzed. The main conclusions are as follows.  (1) Moisture in axisymmetric soil samples migrates from hightemperature region to low-temperature region under the action of the heat source, and its dynamic redistribution causes the change of soil permeability. Accordingly, the reasonable determination of the modified VG model parameters should be paid attention to in the numerical calculation. The analysis shows that the modified VG model parameters of the heating process of four measuring points can take the same value, while the values of the model parameters are different in the cooling process.
(2) The temperature of the heat source directly affects the change of the saturation of the measuring point. The higher the temperature of the heat source, the greater the saturation of the specific position when the soil is thermally stable. However, the change of the soil saturation can be stabilized faster under the lower temperature heat source. Under the same temperature of the heat source, the growth rate of saturation decreases with time, and the bottom saturation of unsaturated soil is greater than the top saturation under gravity. The saturation near the heat source is low. (3) The volumetric strain of unsaturated soil results from the combined action of wet expansion and thermal expansion, which is dominated by moisture expansion near the heat source, and is mainly dominated by thermal expansion far away from the heat source. The volumetric strain near the outer boundary of the soil column is negative, primarily passive deformation, and the volumetric strain in the middle of the soil is positive, primarily active expansion deformation. The change of the total displacement of soil with time is dominated by the z-direction displacement, and its change trend along the radial direction is similar to that in the r-direction. The peak value of the total displacement curve moves outward with time.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
PC; methodology. HZ, GY, ZG, and GY; investigation and research.