Fuzzy Coordination Control Strategy and Thermohydraulic Dynamics Modeling of a Natural Gas Heating System for In Situ Soil Thermal Remediation

Soil contamination remains a global problem. Among the different kinds of remediation technologies, in situ soil thermal remediation has attracted great attention in the environmental field, representing a potential remedial alternative for contaminated soils. Soils need to be heated to a high temperature in thermal remediation, which requires a large amount of energy. For the natural gas heating system in thermal remediation, a fuzzy coordination control strategy and thermohydraulic dynamics model have been proposed in this paper. In order to demonstrate the superiority of the strategy, the other three traditional control strategies are introduced. Analysis of the temperature rise and energy consumption of soils under different control strategies were conducted. The results showed that the energy consumption of fuzzy coordination control strategy is reduced by 33.9% compared to that of the traditional control strategy I, constant natural gas flow and excess air ratio. Further, compared to the traditional control strategy II, constant excess air ratio and desired outlet temperature of wells, the strategy proposed can reduce energy consumption by 48.7%. The results illustrate the superiority of the fuzzy coordination control strategy, and the strategy can greatly reduce energy consumption, thereby reducing the cost of in situ soil thermal remediation.


Introduction
Soil pollution is a major global environmental problem [1,2]. For instance, according to the first national pollution survey in China, released in 2014, 16% of the surveyed land (6.3 million km 2 ) was contaminated beyond acceptable standards [3]. Meanwhile, 342,000 sites of known contamination and a further 2.5 million potentially polluted sites are present in Europe [4]. Due to the severity of soil pollution and the difficulty of its restoration, the remediation of contaminated soil has attracted great attention in the environmental field.
A cornucopia of remediation technologies has been developed to treat contaminated soils [5][6][7][8]. Among these technologies, in situ soil thermal remediation holds an important niche role due to its ability to clean pollutants quickly and reliably [9]. Thermal remediation is a soil remediation process in which heat and vacuum are applied simultaneously to contaminated soils [10]. After the thermal treatment of soils, the vaporized contaminants will be passed through a pollutant treatment system and released into the atmosphere. Then, the treated soil is available for reuse [11]. Studies have indicated that there is no dispersion of pollutants beyond the boundaries of the treated zone in the course of thermal remediation [12]. Meanwhile, TerraTherm Environmental Services, Inc. (TESI) has made a great effort to detect migration in field experiments and demonstrations, and none of these investigations have found evidence of migration after thermal treatment. Therefore, thermal remediation is very effective for treating contaminated soils, and it does not pose a threat to site workers and the community when properly operated. Thermal remediation has been widely applied to the treatment of organic contaminated sites in American and European countries [13][14][15].
Thermal remediation is mainly applied to the removal of volatile, semivolatile organic pollutants, and a small number of volatile inorganic substances such as Hg, As, and Se [16]. Thermal remediation has many advantages over other remediation methods. Firstly, soil temperature can be increased to higher values through thermal remediation, if desired, to values on the order of 600 • C [17]; thus, it can meet cleanup standards for lots of contaminants [18]. Furthermore, thermal remediation is effective for uniformly heating the entire contamination zone [19]. Because thermal conductivity of soils is much less variable, it changes by only a factor of approximately 4 from clay to sand [20]. Further, in situ thermal remediation can be used deep underground or beneath buildings, which would otherwise be difficult or costly to dig up to treat above ground. Thermal remediation can also be used in silty or clayey soils, where other treatment technologies do not perform well. Compared to ex situ, in situ thermal remediation involves no digging and hauling of contaminated soils, generates no dust or odors, minimizes human exposure to hazardous wastes, and is a low-noise operation [21]. Further, in the first full-scale commercial application of in situ thermal remediation, new grass growth reappeared naturally without reseeding after the remediation operation was completed. No damage was observed to the nearby trees [12,22]. Some of the literature has focused on the changes of soil properties after thermal treatment [11,23], and it was found that treated soils have a high abundance and diversity of microorganisms and fauna, but enzymatic activities are weak. Thus, the treated soils may still be suitable for sustaining vegetation, though likely at a slightly diminished capacity when compared with native soils. On the other hand, in situ thermal remediation may be more ecofriendly in sensitive ecosystems because of the lack of soil disturbance [24].
According to the source of energy in thermal remediation, it can be divided into natural gas heating and electric heating. Natural gas heating is more energy efficient than electric heating, but it is accompanied by problems with the control of natural gas flow. A great deal of effort has gone into developing thermal remediation to remove or reduce contaminants in soils. As soil is heated, contaminants in soils are vaporized or destroyed by a number of mechanisms. Heron et al. [25] presented the largest thermal remediation project and the associated challenges of constructing and treating quickly. In this project, an estimated 13,400 kg of chlorinated volatile organic compounds (CVOC) mass was removed, and all soil goals were met. Chang et al. [26] developed the speciation model for mercury over a range of environmental conditions to identify distribution of dissolved mercury species and potential transformations of mercury during thermal remediation. Falciglia et al. [27] used a laboratory scale apparatus to treat five soil size aggregate fractions, and the effect of soil texture on contaminant adsorption and removal was investigated.
Based on the systematic parametric studies mentioned above, soil temperature was recognized as the key parameter for thermal remediation [28]. Because soils have relatively low heat capacity, the initial heating of the site may require long periods of energy input before contaminants vaporize. During the Entropy 2019, 21, 971 3 of 30 thermal remediation process, heating contaminated soils to high temperatures is energy-intensive and, thus, a relatively costly endeavor [24,29]. The treatment costs for soils contaminated with semivolatile organic compounds can be high as $550-770 per metric ton [30]. However, in most cases, the studies focus on the chemical modification of contaminants and the removal of pollutants [31][32][33], and the question of how to reduce energy consumption through an appropriate control strategy is ignored. At present, the heating system of thermal remediation is operated under a simple control strategy that leads to a large amount of energy consumed.
In this paper, a fuzzy coordination control strategy and thermohydraulic dynamics model of natural gas heating systems are proposed for in situ soil thermal remediation. In order to verify the superiority of the strategy, the other three traditional control strategies are introduced. Analysis of the temperature rise and energy consumption of soils under different control strategies was conducted. In Section 2 of this paper, the thermohydraulic models and fuzzy coordination control strategy of heating system are proposed. Section 3 presents the numerical analysis of heating process of soils under four different control strategies. Finally, based on the numerical study results, Section 4 presents the conclusions.

System Description and Dynamical Modeling
In the process of thermal remediation, the entire treatment zone is heated by an array of natural gas heating devices, as Figure 1 shows. The device is composed of a combustor and a thermal well, and the combustor is fixed on the top of wells. The natural gas and air burn in the combustor. Then, the high-temperature flue gases generated flow into the inner pipe of thermal wells. The gases continue to flow downward to the bottom of thermal wells and enter the annular space. Then, the flue gases flow upward away from thermal wells. In thermal remediation, heat is transferred to the outer wall of thermal well's outer pipe by convective heat transfer and heat conduction. After that, soil is heated to high temperatures to achieve the removal of pollutants. The objective of this paper is to propose a control strategy to accomplish energy-saving heating of soils by adjusting the natural gas flow and air flow during the heating process. intensive and, thus, a relatively costly endeavor [24,29]. The treatment costs for soils contaminated with semivolatile organic compounds can be high as $550-770 per metric ton [30]. However, in most cases, the studies focus on the chemical modification of contaminants and the removal of pollutants [31][32][33], and the question of how to reduce energy consumption through an appropriate control strategy is ignored. At present, the heating system of thermal remediation is operated under a simple control strategy that leads to a large amount of energy consumed. In this paper, a fuzzy coordination control strategy and thermohydraulic dynamics model of natural gas heating systems are proposed for in situ soil thermal remediation. In order to verify the superiority of the strategy, the other three traditional control strategies are introduced. Analysis of the temperature rise and energy consumption of soils under different control strategies was conducted. In Section 2 of this paper, the thermohydraulic models and fuzzy coordination control strategy of heating system are proposed. Section 3 presents the numerical analysis of heating process of soils under four different control strategies. Finally, based on the numerical study results, Section 4 presents the conclusions.

System Description and Dynamical Modeling
In the process of thermal remediation, the entire treatment zone is heated by an array of natural gas heating devices, as Figure 1 shows. The device is composed of a combustor and a thermal well, and the combustor is fixed on the top of wells. The natural gas and air burn in the combustor. Then, the high-temperature flue gases generated flow into the inner pipe of thermal wells. The gases continue to flow downward to the bottom of thermal wells and enter the annular space. Then, the flue gases flow upward away from thermal wells. In thermal remediation, heat is transferred to the outer wall of thermal well's outer pipe by convective heat transfer and heat conduction. After that, soil is heated to high temperatures to achieve the removal of pollutants. The objective of this paper is to propose a control strategy to accomplish energy-saving heating of soils by adjusting the natural gas flow and air flow during the heating process.

The Whole Frame of Control Strategy
Considering the heat leakage to surrounding unheated zones and the moisture migration during the heating process in soils, the requirements of energy at different locations are different. At the

The Whole Frame of Control Strategy
Considering the heat leakage to surrounding unheated zones and the moisture migration during the heating process in soils, the requirements of energy at different locations are different. At the same time, the energy demand will be variable when soil temperature changes. On account of these theories, the fuzzy coordination control strategy is proposed. The schematic of control strategy is shown in Figure 2, consisting of a management layer and a control layer.
Entropy 2019, 21, x 4 of 31 same time, the energy demand will be variable when soil temperature changes. On account of these theories, the fuzzy coordination control strategy is proposed. The schematic of control strategy is shown in Figure 2, consisting of a management layer and a control layer. The management layer is a fuzzy system, which puts out the expected outlet temperature of combustors ' fi T and thermal wells ' fo T in accordance with the position of wells w s , soil temperature s T , and water content l w . In this layer, the fuzzy sets, membership functions of variables, operation rules, and fuzzy rules will be determined. The control layer contains multiple heating devices inserted into soils. There are two proportional-integral-derivative (PID) controllers, a combustor and a thermal well in each device. Two PID controllers adjust the natural gas flow and air flow separately. In order to analyze the superiority of control strategy, the mathematical model of heating field, thermal wells, and combustors is developed in this layer.  The management layer is a fuzzy system, which puts out the expected outlet temperature of combustors T f i and thermal wells T f o in accordance with the position of wells s w , soil temperature T s , and water content w l . In this layer, the fuzzy sets, membership functions of variables, operation rules, and fuzzy rules will be determined.
The control layer contains multiple heating devices inserted into soils. There are two proportionalintegral-derivative (PID) controllers, a combustor and a thermal well in each device. Two PID controllers adjust the natural gas flow and air flow separately. In order to analyze the superiority of control strategy, the mathematical model of heating field, thermal wells, and combustors is developed in this layer.

Model of Heating Field
Soil is a typical porous media containing a solid phase, liquid phase, and gas phase. There are pores among soil particles, which are filled with liquid and gas. Studies have shown that the liquid phase and gas phase in soils will migrate from high temperature zone to low temperature zone under the heat drive, which would cause the redistribution of the soil moisture field. The essence of moisture transfer under non-isothermal conditions is the migration of energy, so moisture transfer also affects the change of soil temperature field. Heat and moisture transfer in soils under the temperature gradient and humidity gradient are mainly analyzed in this study.
In thermal remediation, soils are heated by an array of thermal wells. Therefore, the whole site is divided into different unit blocks according to the arrangement of thermal wells, as Figure 3 shows. In the lumped parameter method, it is assumed that soil temperature, water content, and other physical parameters at different locations in a certain unit block are the same, and this method puts emphasis on the difference of heating processes in different blocks. (i, j) indicates the number of each unit block. Soil is a typical porous media containing a solid phase, liquid phase, and gas phase. There are pores among soil particles, which are filled with liquid and gas. Studies have shown that the liquid phase and gas phase in soils will migrate from high temperature zone to low temperature zone under the heat drive, which would cause the redistribution of the soil moisture field. The essence of moisture transfer under non-isothermal conditions is the migration of energy, so moisture transfer also affects the change of soil temperature field. Heat and moisture transfer in soils under the temperature gradient and humidity gradient are mainly analyzed in this study.
In thermal remediation, soils are heated by an array of thermal wells. Therefore, the whole site is divided into different unit blocks according to the arrangement of thermal wells, as Figure 3 shows. In the lumped parameter method, it is assumed that soil temperature, water content, and other physical parameters at different locations in a certain unit block are the same, and this method puts emphasis on the difference of heating processes in different blocks. ( )  In order to simulate the heat and moisture transfer in unsaturated soils, the following simplifying assumptions were made in this study.
(1) Soil is homogeneous and its type does not change along thermal wells, while soil is considered, in a real situation, a non-homogeneous and non-isotropic porous material. The effect of this assumption on heat conduction can be negligible because thermal conductivities of different dry soils are much less variable. By contrast, fluid flow permeabilities of different layers may vary much [12], so the permeabilities of soils at different locations are actually a little variable. (2) Convection of fluid in porous media satisfies Darcy's law. In order to simulate the heat and moisture transfer in unsaturated soils, the following simplifying assumptions were made in this study.
(1) Soil is homogeneous and its type does not change along thermal wells, while soil is considered, in a real situation, a non-homogeneous and non-isotropic porous material. The effect of this assumption on heat conduction can be negligible because thermal conductivities of different dry soils are much less variable. By contrast, fluid flow permeabilities of different layers may vary much [12], so the permeabilities of soils at different locations are actually a little variable. (2) Convection of fluid in porous media satisfies Darcy's law.
(3) The gas phase in soils includes noncondensable gases (such as air) and water vapor. The influence of dry air on heat and moisture migration is neglected. (4) There is no chemical interaction, and the gas is assumed to be ideal gas. (5) It is assumed that the solid, liquid, and gas phases are continuous in unsaturated soil, separately. (6) It is assumed that the migration of liquid and gas do not affect each other. For the mass conservation, the variable quantity of liquid water in a certain unit block is equal to the difference between the amount of migration from the surrounding units and the amount of inside evaporation. As Figure 4 shows, for a single unit block, the mass balance equation of liquid is: where ρ l is the density of liquid water. V T,(i,j) and w l,(i,j) indicate, respectively, total volume and volume liquid moisture content of the unit (i, j),  For the mass conservation, the variable quantity of liquid water in a certain unit block is equal to the difference between the amount of migration from the surrounding units and the amount of inside evaporation. As Figure 4 shows, for a single unit block, the mass balance equation of liquid is: indicate, respectively, total volume and volume liquid moisture content of the unit ( ) J is the mass of liquid migrated from adjacent units in unit time.  According to the Philip and De Vries model [34], the liquid migration equation in unsaturated soils is: where l j indicates liquid migration mass of unit area in unit time. l K and ψ are hydraulic conductivity and soil water potential, respectively. The moisture infiltration on the surface of soils is negligible. In thermal remediation, there are vacuum wells to remove gases in soils, and the extraction of liquid is ignored, so where j l indicates liquid migration mass of unit area in unit time. K l and ψ are hydraulic conductivity and soil water potential, respectively. The moisture infiltration on the surface of soils is negligible. In thermal remediation, there are vacuum wells to remove gases in soils, and the extraction of liquid is ignored, so J l,(i,j),up = 0. In the heating process, there is water infiltration from the unheated zone. These values are set to invariable in this study, as J l,(i,j),e , J l,(i,j),s , J l,(i,j),n , and J l,(i,j),down .
For the lumped parameter method, the liquid fluxes of unit (i, j) with other adjacent units are: where J l,(i−1, j) , J l,(i+1, j) , J l,(i,j−1) , and J l,(i,j+1) indicate water migration mass from unit (i − 1, j), (i + 1, j), (i, j − 1), and (i, j + 1) in unit time. A sd = SL, A u = S 2 , where S is thermal well spacing and L is thermal well depth. D lw and D lT are called isothermal and thermal water diffusivities, respectively.
∂T . The relationship of soil water potential with water content and temperature can be obtained by Gardner and De Vries [35]. ψ = a w l ε −b exp(γ · (T − 273.15)), where a, b and γ are characteristic parameters of soils. ε indicates the porosity of soils.

(b) Vapor Flow Model
Similarly, for the mass conservation, the variable quantity of vapor in a certain unit block is equal to the sum between the amount of migration from surrounding units and the amount of internal evaporation. For a single unit block, the mass balance equation of vapor is: where ρ v is the density of vapor. w v,(i,j) indicates volume gaseous moisture content of the unit (i, j), J v is the mass of vapor migrated from adjacent units in unit time. E v,(i,j) = −E l,(i,j) . Using w l + w v = ε, the left side of Equation (4) can be written as: Comparing Equations (4), (5), and (1), we have: The mechanism of vapor migration in soil is mainly diffusion, which is a transitional diffusion of Fick and Kundsen diffusion. The vapor migration equation is [36]: where j v indicates vapor migration mass of unit area in unit time. D e is the vapor equivalent diffusivity. 1 . D atm and D kn indicate molecular diffusivity and Knudsen diffusivity, respectively. In thermal remediation, there are vacuum wells to extract gases from soils. Extraction is treated as vapor migration on the upper surface of unit block. Therefore, there are two mechanisms of migration on the upper surface: the mass transfer of the upper surface and the extraction of vacuum wells.
According to continuum fluid dynamics theory, the flow velocity of gases, under the extraction of vacuum wells, can be described as the form of a modified Darcy's law: where k is intrinsic permeability of soils. µ gv and p indicate viscosity and pressure of gas phase, respectively. k rg is relative permeability of gas phase, and it can be obtained through the Van Genuchten-Parker empirical formula [37]. k rg S gv = S 1/2 gv 1 − 1 − S gv 1/m 2m , S gv is the saturation of gas phase, and m is an empirical parameter. In addition, there is mass transfer on the surface of soils, which includes convective mass transfer and diffusion mass transfer. Assuming that the convective mass transfer coefficient is β, the vapor mass of convection is: m vec = β[p a − p s ]. According to Fick's law, the mass of diffusion is In summary, the amount of vapor migration on the soil surface is: where A ve and p ve indicate the extraction area and pressure of the vacuum well. p a and p v are gas pressure in the atmosphere and soils. For the lumped parameter method, the vapor fluxes of unit (i, j) with other adjacent units are: where J v,(i−1, j) , J v,(i+1, j) , J v,(i,j−1) , and J v,(i,j+1) indicate vapor migration mass from unit (i − 1, j), ,n , and J v,(i,j),down indicate vapor migration mass from unheated zones. ρ v,n is the vapor density in unheated zones.

(c) Heat Flow Model
For the energy conservation, the enthalpy change in a certain unit block is equal to the sum of heat generated by thermal wells and the net heat flow through the unit. There are solid, liquid, and gas three phases in soils. Therefore, in addition to thermal conduction, the heat flow also includes heat flux caused by liquid migration and vapor migration. As Figure 5 shows, for a single unit block, the energy balance equation is: where M T,(i,j) , c T,(i,j) , and T s,(i,j) indicate total mass, mean specific heat, and temperature in the unit is the heat flowing into soils from the thermal well, which can be calculated by the model of thermal wells. In thermal remediation, the upper surface of soils will be covered with insulation, which is composed of concrete, insulation bricks, gravel layer, etc. Therefore, the thermal resistance of the heat exchange on the soil surface includes the thermal resistance inside soils, the thermal resistance of insulation, and the thermal resistance of heat convection on the upper surface, and solar radiation should also be considered. Thus: Figure 5. Heat flow of a unit block. φ in,(i,j) indicates the heating power of thermal wells: where P (i,j) is the heat flowing into soils from the thermal well, which can be calculated by the model of thermal wells. φ λ,(i−1, j) , φ λ,(i+1, j) , φ λ,(i,j−1) , and φ λ,(i,j+1) indicate the heat flux by thermal conduction from adjacent units. φ λ,(i,j),e , φ λ,(i,j),w , φ λ,(i,j),s , φ λ,(i,j),n , and φ λ,(i,j),down indicate the heat flux by thermal conduction from unheated zones. They can be calculated on the basis of Fourier's law.
In thermal remediation, the upper surface of soils will be covered with insulation, which is composed of concrete, insulation bricks, gravel layer, etc. Therefore, the thermal resistance of the heat exchange on the soil surface includes the thermal resistance inside soils, the thermal resistance of insulation, and the thermal resistance of heat convection on the upper surface, and solar radiation should also be considered. Thus: where , and R 3 = 1 h a . λ ins and δ ins are thermal conductivity and thickness of the insulation. h a is the convective heat transfer coefficient. φ rad indicates the energy radiated by sun on unit area, and α s is the absorption rate of soils.
φ l and φ v indicate the heat flux through liquid migration and vapor migration, respectively.
where H l and H v indicate the enthalpy of liquid water and vapor, is the latent heat of water. J l and J v can be calculated by Equations (3) and (10). φ eva,(i,j) indicates the energy absorbed by evaporation of water in soils. φ eva,(i,j) = γ H 2 O E l,(i,j) .

(d) Three Phases of Soil Warming
During the heating process of thermal remediation, soil temperature can rise up to a very high value (generally around 350 • C [24]). A large number of experimental studies have shown that the temperature history of soils consists of three phases: heat-up phase, boiling phase, and superheating phase [38], as shown in Figure 6.
is the latent heat of water. l J and v J can be calculated by Equations (3) and (10).  During the heating process of thermal remediation, soil temperature can rise up to a very high value (generally around 350 °C [24]). A large number of experimental studies have shown that the temperature history of soils consists of three phases: heat-up phase, boiling phase, and superheating phase [38], as shown in Figure 6. During the heat-up phase, the soil minerals and fluids (mainly water) are heated to the boiling point of water from the initial temperature. There are liquids and gases in the pores of soils. The evaporation of liquid water and liquid migration from unheated zones is negligible in this period.
During the boiling phase, the soil temperature stays at the boiling point until all the pore water has been boiled off. The energy generated by thermal wells is all converted into the latent heat required for evaporation of pore water. The duration of this phase depends on the amount of water to be boiled. When all the pore water has been vaporized, the dry soil can be superheated.
During the superheating phase, the dry soil is heated to the target temperature. The water infiltrated from surrounding units will be evaporated quickly because of the high temperature.
In short, the mechanisms of the three phases are not the same. In each phase, the liquid flow model, vapor flow model, and heat flow model should be slightly modified.

Model of Thermal Well
In order to obtain the heating power of thermal wells, the model of wells is established. Figure  7 shows the structure of thermal wells, which is composed of an inner pipe and an outer pipe, and an annular space is formed between inner pipe and outer pipe. High-temperature flue gases During the heat-up phase, the soil minerals and fluids (mainly water) are heated to the boiling point of water from the initial temperature. There are liquids and gases in the pores of soils. The evaporation of liquid water and liquid migration from unheated zones is negligible in this period.
During the boiling phase, the soil temperature stays at the boiling point until all the pore water has been boiled off. The energy generated by thermal wells is all converted into the latent heat required for evaporation of pore water. The duration of this phase depends on the amount of water to be boiled. When all the pore water has been vaporized, the dry soil can be superheated.
During the superheating phase, the dry soil is heated to the target temperature. The water infiltrated from surrounding units will be evaporated quickly because of the high temperature.
In short, the mechanisms of the three phases are not the same. In each phase, the liquid flow model, vapor flow model, and heat flow model should be slightly modified.

Model of Thermal Well
In order to obtain the heating power of thermal wells, the model of wells is established. Figure 7 shows the structure of thermal wells, which is composed of an inner pipe and an outer pipe, and an annular space is formed between inner pipe and outer pipe. High-temperature flue gases generated by combustors flow in through the inner pipe and then flow out through the outer pipe, thereby heating soils. In order to simplify the analytical process, it is assumed that the heat transfer in the thermal well is quasi-steady. Further, the heat transfer of the fluid in axial direction is ignored, and the heat transfer in the thermal well is considered to be one-dimensional in the radial direction. It is supposed that heat transfer of fluid in radial direction is converted into the temperature variation along the axial direction. Thus, the temperature relations of flue gases in the inner pipe and the annular space along axial direction are respectively established, and an analytical solution of the thermal well's outlet temperature is obtained. The coordinates are set as shown in Figure 7. The origin of coordinates is located at the intersection of the axis of thermal wells and the ground surface. Assumptions: (1) Soil is homogeneous.
(2) The temperature and speed of the same section of fluid in thermal well are the same.
(3) In the model of thermal wells, the thermophysical parameters of soil are invariable, and moisture transfer is ignored. (4) The energy loss at the junction of the inner pipe and the bottom of the outer pipe is not considered.  A microsegment of thermal wells is shown in Figure 8. The heat transfer of the microsegment through the inner pipe is:   A microsegment of thermal wells is shown in Figure 8. The heat transfer of the microsegment through the inner pipe is: where R t1 = 1 . h 1 and h 2 are convective heat transfer coefficients of the inner pipe's inner side and outer side, respectively. r 1 and r 2 are, respectively, the radius of the inner pipe's inner The enthalpy change of flue gases in the inner pipe is: where 1 u , g ρ , and g c indicate the flow rate, density, and specific heat of flue gases in the inner pipe, separately. The thermal balance: 1 1 dQ dH = , then: A heat balance equation is also established for flue gases in the annular space. The heat exchange between flue gases in the annular space with soils is: λ indicate the thermal conductivity of outer pipe and soils, respectively. 5 r is the effective radius of thermal wells. 5 2 S r = , S indicates the thermal well spacing. s T is the soil temperature at 5 r .
The enthalpy change of flue gases in annular space is: where 2 u indicates the flue gases rate in annular space.
The thermal balance of flue gas in annular space is , then: The enthalpy change of flue gases in the inner pipe is: where u 1 , ρ g , and c g indicate the flow rate, density, and specific heat of flue gases in the inner pipe, separately. The thermal balance: dQ 1 = dH 1 , then: A heat balance equation is also established for flue gases in the annular space. The heat exchange between flue gases in the annular space with soils is: . h 3 is the convective heat transfer coefficient of the outer pipe's inner side. r 3 and r 4 are, respectively, the radius of the outer pipe's inner side and outer side. λ 2 and λ 3 indicate the thermal conductivity of outer pipe and soils, respectively. r 5 is the effective radius of thermal wells. r 5 = S 2 , S indicates the thermal well spacing. T s is the soil temperature at r 5 . The enthalpy change of flue gases in annular space is: where u 2 indicates the flue gases rate in annular space. The thermal balance of flue gas in annular space is −dQ 1 + dQ 2 = dH 2 , then: The flue gas density in the thermal well is assumed to be invariable, so the mass flow in the inner pipe is equal to that of annular space.
Boundary conditions are: where T f i is the inlet temperature of thermal wells, which is assumed to be equal to the outlet temperature of combustors. L indicates the depth of inner pipe. Solving Equations (17) and (20) simultaneously, we get: The outlet temperature of the thermal wells is: Then, the heating power of the thermal well can be calculated as: where:

Model of Combustor
In order to calculate the temperature, flow, and specific heat capacity of flue gases at the entrance of thermal wells, the model of a combustor was developed. The schematic diagram of the combustor is shown in Figure 9. Natural gas and air are burned in combustors to generate high-temperature flue gas. The outlet temperature and flow rate of flue gas can be obtained by establishing the mathematical model of chemical combustion in combustors. The natural gas is regarded as methane (CH4) in this study, and air is the combustion-supporting gas, which is composed of 78% nitrogen, 21% oxygen, and 1% other gases. Further, more air will be provided for the complete combustion of natural gas. In this paper, excess air ratio α is defined as the ratio of excess air to theoretical air. Then, the chemical equation can be obtained: where Q indicates the heat generated by combustion. The flue gas produced by methane combustion is a mixture of different gases. Its thermophysical parameters should be analyzed before the calculation of the flue gas's temperature and flow. All gases are considered as ideal gases during the calculation.
The methane combustion is supposed to be complete. The proportion of each gas The average specific heat of flue gas is: The specific heat capacity of gas is different at different temperatures. The molar specific heat of components in flue gases The density of gas at different temperatures is: According to , the density of flue gas is: where 2 T indicates the outlet Kelvin temperature of combustors.
Input energy is equal to the energy of output according to conservation of energy, that is, the heat of methane combustion is equal to the energy required for the heating up of flue gases and the latent heat demanded by water vaporization. Moreover, there will be a certain loss of heat in combustion, so combustion efficiency ξ is introduced.
where 1 G indicates molar flow of natural gas. 2 G indicates molar flow of flue gases.  The flue gas produced by methane combustion is a mixture of different gases. Its thermophysical parameters should be analyzed before the calculation of the flue gas's temperature and flow. All gases are considered as ideal gases during the calculation.
The methane combustion is supposed to be complete. The proportion of each gas w CO 2 , w H 2 O , w N 2 , and w O 2 can be obtained based on the chemical equation.
The average specific heat of flue gas is: The specific heat capacity of gas is different at different temperatures. The molar specific heat of components in flue gases c CO 2 , c H 2 O , c N 2 , and c O 2 can be obtained according to the literature. Then: . The density of gas at different temperatures is: According to ρ g = ρ CO 2 · w CO 2 + ρ H 2 O · w H 2 O + ρ N 2 · w N 2 + ρ O 2 · w O 2 , the density of flue gas is: where T 2 indicates the outlet Kelvin temperature of combustors. Input energy is equal to the energy of output according to conservation of energy, that is, the heat of methane combustion is equal to the energy required for the heating up of flue gases and the latent heat demanded by water vaporization. Moreover, there will be a certain loss of heat in combustion, so combustion efficiency ξ is introduced.
where G 1 indicates molar flow of natural gas. After Equation (31) is solved, it can be obtained that: where t 1 is the initial temperature of gas, which is usually low, so it is supposed to be 0. Then, the outlet temperature of combustors is: According to Equation (26), the flue gas flow is: where G g is the volume flow of flue gases. G N indicates the natural gas volume flow. T 1 and T 2 are the Kelvin temperature of natural gas and flue gas, separately.

Fuzzy Coordination Control Strategy
The graphical illustration of the fuzzy coordination control strategy is presented in Figure 10. It accepts crisp mathematical values as input, and these crisp values are then converted into fuzzy sets by fuzzification interface [39]. The inference engine uses fuzzy rules presented in the knowledge base and generates output fuzzy sets. These fuzzy sets are then converted back to crisp values by the defuzzification interface. The fuzzy rules are defined by the system developer and vary from problem to problem [40]. where 1 t is the initial temperature of gas, which is usually low, so it is supposed to be 0. Then, the outlet temperature of combustors is: According to Equation (26), the flue gas flow is:

Fuzzy Coordination Control Strategy
The graphical illustration of the fuzzy coordination control strategy is presented in Figure 10. It accepts crisp mathematical values as input, and these crisp values are then converted into fuzzy sets by fuzzification interface [39]. The inference engine uses fuzzy rules presented in the knowledge base and generates output fuzzy sets. These fuzzy sets are then converted back to crisp values by the defuzzification interface. The fuzzy rules are defined by the system developer and vary from problem to problem [40]. The energy requirements of soils are variational at different locations. Thus, the position of thermal well s w is considered to reduce the impact of boundaries. s w indicates the sum of lateral and longitudinal distances from the thermal wells to the well on the border. Soil temperature T s and water content w l are considered to judge the warming phase of soils. The membership functions of input and output variables are given in Figures 11 and 12. The Gaussian function has been used to ensure a smoother curve of output variables.
The output is determined using the fuzzy rules in the following form: IF s w is E i and T s is CE j and w l is CU k , where E i , CE j , CU k , CT (i,j,k) , and CP (i,j,k) are the fuzzy values of s w , T s , w l , T f i , and T f o . There are 27 fuzzy rules in our fuzzy system. Because the water content changes during the boiling phase, the value of water content is considered only at this phase. Table 1 shows the complete list of fuzzy rules, which suggests the following: (1) If the position of wells is closer to the boundary, the set values of combustors' and wells' outlet temperature should be higher. (2) As the soil temperature increases and the water content decreases gradually, the desired outlet temperatures of combustors and thermal wells are set to be larger. The temperatures are set to be lower in the early stage and can reduce energy consumption. There are 27 fuzzy rules in our fuzzy system. Because the water content changes during the boiling phase, the value of water content is considered only at this phase. Table 1 shows the complete list of fuzzy rules, which suggests the following: (1) If the position of wells is closer to the boundary, the set values of combustors' and wells' outlet temperature should be higher. (2) As the soil temperature increases and the water content decreases gradually, the desired outlet temperatures of combustors and thermal wells are set to be larger. The temperatures are set to be lower in the early stage and can reduce energy consumption. There are 27 fuzzy rules in our fuzzy system. Because the water content changes during the boiling phase, the value of water content is considered only at this phase. Table 1 shows the complete list of fuzzy rules, which suggests the following: (1) If the position of wells is closer to the boundary, the set values of combustors' and wells' outlet temperature should be higher. (2) As the soil temperature increases and the water content decreases gradually, the desired outlet temperatures of combustors and thermal wells are set to be larger. The temperatures are set to be lower in the early stage and can reduce energy consumption.

Results and Discussion
The focus of this paper was to propose a coordination control strategy of multiple wells for the treatment zones, so as to reduce the energy consumption during the heating process without changing the duration. In addition to the fuzzy coordination control strategy, the energy consumptions of the field under the other three control strategies were also analyzed.
The mathematical models presented above were built in MATLAB/Simulink, and the calculation results under different control strategies were obtained. Figure 13

Results and Discussion
The focus of this paper was to propose a coordination control strategy of multiple wells for the treatment zones, so as to reduce the energy consumption during the heating process without changing the duration. In addition to the fuzzy coordination control strategy, the energy consumptions of the field under the other three control strategies were also analyzed.
The mathematical models presented above were built in MATLAB/Simulink, and the calculation results under different control strategies were obtained. Figure 13  In this work, energy consumption is defined as the amount of natural gas consumed by the heating of soils. The energy consumption of a single thermal well can be calculated by: where , i j G indicates the energy consumption of a single thermal well. z τ is the heating time.  In this work, energy consumption is defined as the amount of natural gas consumed by the heating of soils. The energy consumption of a single thermal well can be calculated by: where G i,j indicates the energy consumption of a single thermal well. τ z is the heating time. G N,(i,j) is the natural gas flow at different time. Then, the energy consumption of all thermal wells in the entire treatment zone is: where i max and j max are the number of rows and columns of thermal wells' pattern. In this paper, four cases were simulated in the MATLAB software. The structure diagrams of the four control strategies are shown in Figure 14. The differences are: In this work, energy consumption is defined as the amount of natural gas consumed by the heating of soils. The energy consumption of a single thermal well can be calculated by: where , i j G indicates the energy consumption of a single thermal well. z τ is the heating time.  In addition to the control strategy, there are no differences in mathematical models and simulation parameters of four cases. The soil type in simulation is sandy soil; its properties are shown in Table 2.

J/ kg K ⋅
Other simulation parameters are shown in Table 3. It can be seen that 101.25 m 3 soils were treated during the process, which are about 194 tons.

Case 1: Constant Natural Gas Flow and Excess Air Ratio
In case 1, natural gas flow and excess air ratio are constant during the heating process. The values are set to 3 3 2.05 10 m /s − × and 1.5, respectively. Because the soil is assumed to be isotropic, parameters such as temperature and water content of the site are symmetrically distributed, so the units (1,1), (1,2), and (2,2) are mainly concerned. Figure 15 shows the variation of temperatures and In addition to the control strategy, there are no differences in mathematical models and simulation parameters of four cases. The soil type in simulation is sandy soil; its properties are shown in Table 2. Other simulation parameters are shown in Table 3. It can be seen that 101.25 m 3 soils were treated during the process, which are about 194 tons.

Case 1: Constant Natural Gas Flow and Excess Air Ratio
In case 1, natural gas flow and excess air ratio are constant during the heating process. The values are set to 2.05 × 10 −3 m 3 /s and 1.5, respectively. Because the soil is assumed to be isotropic, parameters such as temperature and water content of the site are symmetrically distributed, so the units (1,1), (1,2), and (2,2) are mainly concerned. Figure 15 shows the variation of temperatures and water contents in units (1,1), (1,2), and (2,2). The temperature rise can be divided into three phases: heat-up phase, boiling phase, and superheating phase. The water content declines to 0 during the boiling phase.
The temperature rise curves of three units are quite different. Due to the influence of boundary heat leakage, unit (2,2) reaches the target temperature firstly, unit (1,2) is the second, and unit (1,1) is the last. The duration represents the heating time required to reach the target temperature. The durations of the three units are 58.62, 39.83, and 33.37 d, separately. The heating time of the site is supposed to be the largest duration of all units, and the temperature of unit (2,2) will continue to rise to 800 K, eventually.
water contents in units (1,1), (1,2), and (2,2). The temperature rise can be divided into three phases: heat-up phase, boiling phase, and superheating phase. The water content declines to 0 during the boiling phase.
The temperature rise curves of three units are quite different. Due to the influence of boundary heat leakage, unit (2,2) reaches the target temperature firstly, unit (1,2) is the second, and unit (1,1) is the last. The duration represents the heating time required to reach the target temperature. The durations of the three units are 58.62, 39.83, and 33.37 d, separately. The heating time of the site is supposed to be the largest duration of all units, and the temperature of unit (2,2) will continue to rise to 800 K, eventually.  Figure 16 illustrates the variation of natural gas flow, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). In order to analyze the changes of parameters perfectly during the heating, there is no natural gas burning in the first half day of the simulation. In this case, the natural gas flow and excess air ratio of all units are constant, as Figure  16a shows. The outlet temperature of combustors is only related to the excess air ratio. Thus, it is invariable since the excess air ratio is constant, and it can be seen that it is about 1080 K from Figure  16c. Figure 16b represents that the heating power of thermal wells decreases over time. From Figure  16d, it can be observed that the outlet temperature of thermal wells increases over time. Soil thermal conductivity decreases as soil moisture content reduces. As the soil temperature increases and thermal conductivity decreases, the amount of heat absorbed by soils becomes less gradual. Therefore, when natural gas flow is constant, the outlet temperature of thermal wells will gradually increase. Meanwhile, because the soil temperature of unit (2, 2) is higher, there is a higher outlet temperature in its thermal well, and its heating power is lower.  Figure 16 illustrates the variation of natural gas flow, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). In order to analyze the changes of parameters perfectly during the heating, there is no natural gas burning in the first half day of the simulation. In this case, the natural gas flow and excess air ratio of all units are constant, as Figure 16a shows. The outlet temperature of combustors is only related to the excess air ratio. Thus, it is invariable since the excess air ratio is constant, and it can be seen that it is about 1080 K from Figure 16c. Figure 16b represents that the heating power of thermal wells decreases over time. From Figure 16d, it can be observed that the outlet temperature of thermal wells increases over time. Soil thermal conductivity decreases as soil moisture content reduces. As the soil temperature increases and thermal conductivity decreases, the amount of heat absorbed by soils becomes less gradual. Therefore, when natural gas flow is constant, the outlet temperature of thermal wells will gradually increase. Meanwhile, because the soil temperature of unit (2, 2) is higher, there is a higher outlet temperature in its thermal well, and its heating power is lower.
heat-up phase, boiling phase, and superheating phase. The water content declines to 0 during the boiling phase.
The temperature rise curves of three units are quite different. Due to the influence of boundary heat leakage, unit (2,2) reaches the target temperature firstly, unit (1,2) is the second, and unit (1,1) is the last. The duration represents the heating time required to reach the target temperature. The durations of the three units are 58.62, 39.83, and 33.37 d, separately. The heating time of the site is supposed to be the largest duration of all units, and the temperature of unit (2,2) will continue to rise to 800 K, eventually.  Figure 16 illustrates the variation of natural gas flow, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). In order to analyze the changes of parameters perfectly during the heating, there is no natural gas burning in the first half day of the simulation. In this case, the natural gas flow and excess air ratio of all units are constant, as Figure  16a shows. The outlet temperature of combustors is only related to the excess air ratio. Thus, it is invariable since the excess air ratio is constant, and it can be seen that it is about 1080 K from Figure  16c. Figure 16b represents that the heating power of thermal wells decreases over time. From Figure  16d, it can be observed that the outlet temperature of thermal wells increases over time. Soil thermal conductivity decreases as soil moisture content reduces. As the soil temperature increases and thermal conductivity decreases, the amount of heat absorbed by soils becomes less gradual. Therefore, when natural gas flow is constant, the outlet temperature of thermal wells will gradually increase. Meanwhile, because the soil temperature of unit (2, 2) is higher, there is a higher outlet temperature in its thermal well, and its heating power is lower. The energy consumption of entire treatment zones is 92673 m 3 in case 1. The proportions of energy consumption during each phase are shown in Figure 17. It can be seen that the superheating phase's energy consumption is largest, and the heat-up phase's energy consumption is very little. In this case, the natural gas flow remains unchanged. Thus, the proportion of energy consumption in different phases is mainly related to the duration of each phase.

Case 2: Constant Excess Air Ratio and Desired Outlet Temperature of Thermal Wells
In case 2, the excess air ratio and desired outlet temperature of thermal wells are constant during the heating process. They are set to 1.5 and 1000 K, respectively. The outlet temperature of thermal wells is taken as the controlled variable, and 1000 K is the set value. Further, the natural gas flow is adjusted by the deviation between the true value of the controlled variable and set value. This control method is most commonly used in engineering. Figure 18 illustrates the variation of soil temperature and water content in case 2. The curves are very similar to those of case 1; the soil temperature rise can also be divided into three phases. The durations of units (1,1), (1,2), and (2,2) are 58.25, 36.05, and 29.94 d, separately. As soil temperature increases, the heating power decreases gradually. At the same time, the heat leakage also increases when temperature rises. Therefore, the warming trend of soils will gradually slow down in the superheating phase, which is especially obvious in unit (1,1) because of the effect of unheated zones. The energy consumption of entire treatment zones is 92673 m 3 in case 1. The proportions of energy consumption during each phase are shown in Figure 17. It can be seen that the superheating phase's energy consumption is largest, and the heat-up phase's energy consumption is very little. In this case, the natural gas flow remains unchanged. Thus, the proportion of energy consumption in different phases is mainly related to the duration of each phase. The energy consumption of entire treatment zones is 92673 m 3 in case 1. The proportions of energy consumption during each phase are shown in Figure 17. It can be seen that the superheating phase's energy consumption is largest, and the heat-up phase's energy consumption is very little. In this case, the natural gas flow remains unchanged. Thus, the proportion of energy consumption in different phases is mainly related to the duration of each phase.

Case 2: Constant Excess Air Ratio and Desired Outlet Temperature of Thermal Wells
In case 2, the excess air ratio and desired outlet temperature of thermal wells are constant during the heating process. They are set to 1.5 and 1000 K, respectively. The outlet temperature of thermal wells is taken as the controlled variable, and 1000 K is the set value. Further, the natural gas flow is adjusted by the deviation between the true value of the controlled variable and set value. This control method is most commonly used in engineering. Figure 18 illustrates the variation of soil temperature and water content in case 2. The curves are very similar to those of case 1; the soil temperature rise can also be divided into three phases. The durations of units (1,1), (1,2), and (2,2) are 58.25, 36.05, and 29.94 d, separately. As soil temperature increases, the heating power decreases gradually. At the same time, the heat leakage also increases when temperature rises. Therefore, the warming trend of soils will gradually slow down in the superheating phase, which is especially obvious in unit (1,1) because of the effect of unheated zones.

Case 2: Constant Excess Air Ratio and Desired Outlet Temperature of Thermal Wells
In case 2, the excess air ratio and desired outlet temperature of thermal wells are constant during the heating process. They are set to 1.5 and 1000 K, respectively. The outlet temperature of thermal wells is taken as the controlled variable, and 1000 K is the set value. Further, the natural gas flow is adjusted by the deviation between the true value of the controlled variable and set value. This control method is most commonly used in engineering. Figure 18 illustrates the variation of soil temperature and water content in case 2. The curves are very similar to those of case 1; the soil temperature rise can also be divided into three phases. The durations of units (1,1), (1,2), and (2,2) are 58.25, 36.05, and 29.94 d, separately. As soil temperature increases, the heating power decreases gradually. At the same time, the heat leakage also increases when temperature rises. Therefore, the warming trend of soils will gradually slow down in the superheating phase, which is especially obvious in unit (1,1) because of the effect of unheated zones.  Figure 19 represents the variation of natural gas flow, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). In this case, the excess air ratio and desired outlet temperature of thermal wells are set to be unchanged, which can be seen in Figure  19c,d. As the soil temperature rises and thermal conductivity decreases, the amount of heat absorbed by soils becomes less gradual. Therefore, when outlet temperature of thermal wells is constant, the natural gas flow will gradually decrease, as Figure 19a shows. Meanwhile, because the soil temperature of unit (2,2) is highest, the natural gas flow is littlest in unit (2,2). From Figure 19b, it can be seen that the heating power decreases over time, and the value of unit (2,2) is the smallest.   Figure 19 represents the variation of natural gas flow, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). In this case, the excess air ratio and desired outlet temperature of thermal wells are set to be unchanged, which can be seen in Figure 19c,d. As the soil temperature rises and thermal conductivity decreases, the amount of heat absorbed by soils becomes less gradual. Therefore, when outlet temperature of thermal wells is constant, the natural gas flow will gradually decrease, as Figure 19a shows. Meanwhile, because the soil temperature of unit (2,2) is highest, the natural gas flow is littlest in unit (2,2). From Figure 19b, it can be seen that the heating power decreases over time, and the value of unit (2,2) is the smallest.  Figure 19 represents the variation of natural gas flow, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). In this case, the excess air ratio and desired outlet temperature of thermal wells are set to be unchanged, which can be seen in Figure  19c,d. As the soil temperature rises and thermal conductivity decreases, the amount of heat absorbed by soils becomes less gradual. Therefore, when outlet temperature of thermal wells is constant, the natural gas flow will gradually decrease, as Figure 19a shows. Meanwhile, because the soil temperature of unit (2,2) is highest, the natural gas flow is littlest in unit (2,2). From Figure 19b, it can be seen that the heating power decreases over time, and the value of unit (2,2) is the smallest.  the boiling phase and superheating phase. The energy consumption of the boiling phase can account for nearly 40%, so the moisture content has a great impact on the temperature rise of soils. Therefore, it will be necessary to control the water influx with temporary bulkheads, freeze walls, or de-watering in some sites. The energy consumption of entire treatment zones is 119250 m 3 in case 2. The proportions of each phases are shown in Figure 20. In the process of thermal remediation, a mass of energy is consumed in the boiling phase and superheating phase. The energy consumption of the boiling phase can account for nearly 40%, so the moisture content has a great impact on the temperature rise of soils. Therefore, it will be necessary to control the water influx with temporary bulkheads, freeze walls, or de-watering in some sites.

Case 3: Constant Desired Outlet Temperature of Combustors and Thermal Wells
In case 3, the desired outlet temperatures of combustors and thermal wells are constant. The values are set to be 1100 and 1004 K, respectively. The flow of natural gas and air is adjusted according to the deviation of true value and set value of two temperatures, respectively. Figure 21 illustrates the variation of soil temperature and water content in case 3. The curves are the same as those of case 1 and case 2. The durations of units (1,1), (1,2), and (2,2) are 58.43, 36.14, and 29.99 d.  Figure 22 represents the variation of natural gas flow, excess air ratio, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). In this case, the desired outlet temperatures of combustors and thermal wells are set to be unchanged, which can be seen in Figure 22d,e. The outlet temperature of combustors is constant, so the excess air ratio is invariable over time, as Figure 22b shows. From Figure 22a,c, it can be seen that the natural gas flow and heating power gradually decreases over time, like that of case 2.

Case 3: Constant Desired Outlet Temperature of Combustors and Thermal Wells
In case 3, the desired outlet temperatures of combustors and thermal wells are constant. The values are set to be 1100 and 1004 K, respectively. The flow of natural gas and air is adjusted according to the deviation of true value and set value of two temperatures, respectively. Figure 21 illustrates the variation of soil temperature and water content in case 3. The curves are the same as those of case 1 and case 2. The durations of units (1,1), (1,2), and (2,2) are 58.43, 36.14, and 29.99 d. The energy consumption of entire treatment zones is 119250 m 3 in case 2. The proportions of each phases are shown in Figure 20. In the process of thermal remediation, a mass of energy is consumed in the boiling phase and superheating phase. The energy consumption of the boiling phase can account for nearly 40%, so the moisture content has a great impact on the temperature rise of soils. Therefore, it will be necessary to control the water influx with temporary bulkheads, freeze walls, or de-watering in some sites.

Case 3: Constant Desired Outlet Temperature of Combustors and Thermal Wells
In case 3, the desired outlet temperatures of combustors and thermal wells are constant. The values are set to be 1100 and 1004 K, respectively. The flow of natural gas and air is adjusted according to the deviation of true value and set value of two temperatures, respectively. Figure 21 illustrates the variation of soil temperature and water content in case 3. The curves are the same as those of case 1 and case 2. The durations of units (1,1), (1,2), and (2,2) are 58.43, 36.14, and 29.99 d.  Figure 22 represents the variation of natural gas flow, excess air ratio, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). In this case, the desired outlet temperatures of combustors and thermal wells are set to be unchanged, which can be seen in Figure 22d,e. The outlet temperature of combustors is constant, so the excess air ratio is invariable over time, as Figure 22b shows. From Figure 22a,c, it can be seen that the natural gas flow and heating power gradually decreases over time, like that of case 2.  Figure 22 represents the variation of natural gas flow, excess air ratio, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). In this case, the desired outlet temperatures of combustors and thermal wells are set to be unchanged, which can be seen in Figure 22d,e. The outlet temperature of combustors is constant, so the excess air ratio is invariable over time, as Figure 22b shows. From Figure 22a,c, it can be seen that the natural gas flow and heating power gradually decreases over time, like that of case 2.

Case 4: Fuzzy Coordination Control Strategy
In case 4, the heating process of soils is controlled by the fuzzy coordination control strategy. The desired outlet temperature of combustors and thermal wells are adjusted by the fuzzy system. Figure 24 illustrates the variation of soil temperature and water content in case 4. It can be seen that the difference among three units' temperature rise curves is very little after the adjustment of fuzzy system, and the three units reach the target temperature almost simultaneously. The durations of units (1,1), (1,2), and (2,2) are 58.75, 58.37, and 56.94 d.

Case 4: Fuzzy Coordination Control Strategy
In case 4, the heating process of soils is controlled by the fuzzy coordination control strategy. The desired outlet temperature of combustors and thermal wells are adjusted by the fuzzy system. Figure 24 illustrates the variation of soil temperature and water content in case 4. It can be seen that the difference among three units' temperature rise curves is very little after the adjustment of fuzzy system, and the three units reach the target temperature almost simultaneously. The durations of units (1,1), (1,2), and (2,2) are 58.75, 58.37, and 56.94 d.  Figure 25 represents the variation of natural gas flow, excess air ratio, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). From Figure  25d,e, it can be observed that the outlet temperatures of combustors and thermal wells increase gradually with time after the adjustment of fuzzy system, and the closer the distance is to the boundary, the larger the values are. Therefore, the set values of unit (1,1) are largest among three units. The natural gas flow increases over time, and the excess air ratio reduces over time, as Figure  25a,b shows. From Figure 25c, it can be seen that the heating power keeps a small fluctuation. Since unit (1,1) has the largest amount of heat dissipation, its heating power is the highest.  Figure 25 represents the variation of natural gas flow, excess air ratio, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). From Figure 25d,e, it can be observed that the outlet temperatures of combustors and thermal wells increase gradually with time after the adjustment of fuzzy system, and the closer the distance is to the boundary, the larger the values are. Therefore, the set values of unit (1,1) are largest among three units. The natural gas flow increases over time, and the excess air ratio reduces over time, as Figure 25a,b shows. From Figure 25c, it can be seen that the heating power keeps a small fluctuation. Since unit (1,1) has the largest amount of heat dissipation, its heating power is the highest.

Case 4: Fuzzy Coordination Control Strategy
In case 4, the heating process of soils is controlled by the fuzzy coordination control strategy. The desired outlet temperature of combustors and thermal wells are adjusted by the fuzzy system. Figure 24 illustrates the variation of soil temperature and water content in case 4. It can be seen that the difference among three units' temperature rise curves is very little after the adjustment of fuzzy system, and the three units reach the target temperature almost simultaneously. The durations of units (1,1), (1,2), and (2,2) are 58.75, 58.37, and 56.94 d.  Figure 25 represents the variation of natural gas flow, excess air ratio, heating power, and outlet temperature of combustors and thermal wells over time for units (1,1), (1,2), and (2,2). From Figure  25d,e, it can be observed that the outlet temperatures of combustors and thermal wells increase gradually with time after the adjustment of fuzzy system, and the closer the distance is to the boundary, the larger the values are. Therefore, the set values of unit (1,1) are largest among three units. The natural gas flow increases over time, and the excess air ratio reduces over time, as Figure  25a,b shows. From Figure 25c, it can be seen that the heating power keeps a small fluctuation. Since unit (1,1) has the largest amount of heat dissipation, its heating power is the highest.

Comparison of Energy Consumption in Four Cases
After the simulation analysis of four cases, the energy consumption of different cases can be compared and analyzed to verify the superiority of the proposed control strategy for reducing energy consumption. Figure 27 illustrates the total energy consumption of the entire treatment zones in the four cases. It can be seen that the energy consumption of case 2 is the highest, 119250 m 3 , and the energy consumption of case 4 is the lowest, 61217 m 3 . Compared to that of case 1, the total energy consumption of cases 2 and 3 increased by 28.7% and 14.1%, respectively, while that of case 4 decreased by 33.9%. Meanwhile, the durations of the four cases are not much different, all being around 58 days. According to research, the control method of case 2 is commonly used by some companies, and compared to this method, the energy consumption of case 4 can be decreased by 48.7%. Therefore, the fuzzy coordination control strategy proposed in this paper can greatly reduce energy consumption without changing the total duration and achieve the energy-saving heating of soils.

Comparison of Energy Consumption in Four Cases
After the simulation analysis of four cases, the energy consumption of different cases can be compared and analyzed to verify the superiority of the proposed control strategy for reducing energy consumption. Figure 27 illustrates the total energy consumption of the entire treatment zones in the four cases. It can be seen that the energy consumption of case 2 is the highest, 119250 m 3 , and the energy consumption of case 4 is the lowest, 61217 m 3 . Compared to that of case 1, the total energy consumption of cases 2 and 3 increased by 28.7% and 14.1%, respectively, while that of case 4 decreased by 33.9%. Meanwhile, the durations of the four cases are not much different, all being around 58 days. According to research, the control method of case 2 is commonly used by some companies, and compared to this method, the energy consumption of case 4 can be decreased by 48.7%. Therefore, the fuzzy coordination control strategy proposed in this paper can greatly reduce energy consumption without changing the total duration and achieve the energy-saving heating of soils.

Comparison of Energy Consumption in Four Cases
After the simulation analysis of four cases, the energy consumption of different cases can be compared and analyzed to verify the superiority of the proposed control strategy for reducing energy consumption. Figure 27 illustrates the total energy consumption of the entire treatment zones in the four cases. It can be seen that the energy consumption of case 2 is the highest, 119,250 m 3 , and the energy consumption of case 4 is the lowest, 61,217 m 3 . Compared to that of case 1, the total energy consumption of cases 2 and 3 increased by 28.7% and 14.1%, respectively, while that of case 4 decreased by 33.9%. Meanwhile, the durations of the four cases are not much different, all being around 58 days. According to research, the control method of case 2 is commonly used by some companies, and compared to this method, the energy consumption of case 4 can be decreased by 48.7%. Therefore, the fuzzy coordination control strategy proposed in this paper can greatly reduce energy consumption without changing the total duration and achieve the energy-saving heating of soils. Entropy 2019, 21, x 28 of 31 Figure 27. The total energy consumption of four cases. Figure 28 represents the energy consumption of each unit in four cases. Compared to that in case 1, the energy consumption of units (1,1), (1,2), and (2,2) in case 4 reduced by 13.2%, 46.9%, and 65.0%, separately. Further, among the three units, the degree of reduction in unit (2,2) is the highest. This is because in cases 1, 2, and 3, the difference in temperature rise of different units is large, and unit (2,2) will reach the target temperature as soon as possible. However, the entire site has not reached the target temperature, and the unit is still heated, so that its temperature will continue to rise, close to 800 K. Further, because of the large temperature difference among units, a considerable part of the energy will be transmitted to other units from unit (2,2), so the energy consumption is relatively high. In case 4, the temperature difference among different units is small, and the durations of different units are not much different, with less heat removed from unit (2,2). Therefore, the reduction proportion of unit (2,2) is the largest, and that of unit (1,1) is the littlest.

Conclusions
In order to reduce the energy consumption of in situ soil thermal remediation, this paper presented a thermohydraulic dynamics model and fuzzy coordination control strategy for a natural gas heating system. Multiple heating systems are controlled coordinately to heat 101.25 m 3 soils in the treatment site. In addition to the fuzzy coordination control strategy, the other three traditional control strategies were also analyzed and compared. Several conclusions can be obtained: (1) When the site is heated by indiscriminate energy, the soil at the boundary is warmed more slowly due to the larger heat dissipation; (2) In the three phases of soil warming, the proportions of energy consumption in the boiling phase and superheating phase are the largest, and that of the heat-up phase is the smallest; (3) Compared to the traditional control strategy I, that is, constant natural gas flow and excess air ratio, the fuzzy coordination control strategy can reduce energy consumption by 33.9%. (4) Compared to the traditional control strategy II, that is, constant excess air ratio and desired outlet temperature of wells, the fuzzy coordination control strategy can reduce energy consumption by 48.7%.  Figure 28 represents the energy consumption of each unit in four cases. Compared to that in case 1, the energy consumption of units (1,1), (1,2), and (2,2) in case 4 reduced by 13.2%, 46.9%, and 65.0%, separately. Further, among the three units, the degree of reduction in unit (2,2) is the highest. This is because in cases 1, 2, and 3, the difference in temperature rise of different units is large, and unit (2,2) will reach the target temperature as soon as possible. However, the entire site has not reached the target temperature, and the unit is still heated, so that its temperature will continue to rise, close to 800 K. Further, because of the large temperature difference among units, a considerable part of the energy will be transmitted to other units from unit (2,2), so the energy consumption is relatively high. In case 4, the temperature difference among different units is small, and the durations of different units are not much different, with less heat removed from unit (2,2). Therefore, the reduction proportion of unit (2,2) is the largest, and that of unit (1,1) is the littlest.  Figure 28 represents the energy consumption of each unit in four cases. Compared to that in case 1, the energy consumption of units (1,1), (1,2), and (2,2) in case 4 reduced by 13.2%, 46.9%, and 65.0%, separately. Further, among the three units, the degree of reduction in unit (2,2) is the highest. This is because in cases 1, 2, and 3, the difference in temperature rise of different units is large, and unit (2,2) will reach the target temperature as soon as possible. However, the entire site has not reached the target temperature, and the unit is still heated, so that its temperature will continue to rise, close to 800 K. Further, because of the large temperature difference among units, a considerable part of the energy will be transmitted to other units from unit (2,2), so the energy consumption is relatively high. In case 4, the temperature difference among different units is small, and the durations of different units are not much different, with less heat removed from unit (2,2). Therefore, the reduction proportion of unit (2,2) is the largest, and that of unit (1,1) is the littlest.

Conclusions
In order to reduce the energy consumption of in situ soil thermal remediation, this paper presented a thermohydraulic dynamics model and fuzzy coordination control strategy for a natural gas heating system. Multiple heating systems are controlled coordinately to heat 101.25 m 3 soils in the treatment site. In addition to the fuzzy coordination control strategy, the other three traditional control strategies were also analyzed and compared. Several conclusions can be obtained: (1) When the site is heated by indiscriminate energy, the soil at the boundary is warmed more slowly due to the larger heat dissipation; (2) In the three phases of soil warming, the proportions of energy consumption in the boiling phase and superheating phase are the largest, and that of the heat-up phase is the smallest; (3) Compared to the traditional control strategy I, that is, constant natural gas flow and excess air ratio, the fuzzy coordination control strategy can reduce energy consumption by 33.9%. (4) Compared to the traditional control strategy II, that is, constant excess air ratio and desired outlet temperature of wells, the fuzzy coordination control strategy can reduce energy consumption by 48.7%.

Conclusions
In order to reduce the energy consumption of in situ soil thermal remediation, this paper presented a thermohydraulic dynamics model and fuzzy coordination control strategy for a natural gas heating system. Multiple heating systems are controlled coordinately to heat 101.25 m 3 soils in the treatment site. In addition to the fuzzy coordination control strategy, the other three traditional control strategies were also analyzed and compared. Several conclusions can be obtained: (1) When the site is heated by indiscriminate energy, the soil at the boundary is warmed more slowly due to the larger heat dissipation; (2) In the three phases of soil warming, the proportions of energy consumption in the boiling phase and superheating phase are the largest, and that of the heat-up phase is the smallest; (3) Compared to the traditional control strategy I, that is, constant natural gas flow and excess air ratio, the fuzzy coordination control strategy can reduce energy consumption by 33.9%.
(4) Compared to the traditional control strategy II, that is, constant excess air ratio and desired outlet temperature of wells, the fuzzy coordination control strategy can reduce energy consumption by 48.7%. (5) After the regulation of the fuzzy coordination control strategy, the soil heating rate is almost the same at different locations, and compared to soils in different locations, the energy consumption of soils located centrally was mostly reduced.
The results showed that the fuzzy coordination control strategy can greatly reduce energy consumption, thereby reducing the cost of in situ soil thermal remediation. This is conducive to the promotion and application of thermal remediation for the treatment of contaminated soils. Further, the models and methodology are valuable for different types of soils and different project requirements, which has a very important guiding significance for engineering application of thermal remediation.