Planned Heating Control Strategy and Thermodynamic Modeling of a Natural Gas Thermal Desorption System for Contaminated Soil

This paper presents a planned heating control strategy applied for a natural gas thermal desorption system for polluted soil to achieve the dynamic adjustment of the heating time and energy consumption. A lumped-parameter model for the proposed system is established to examine effects of the natural gas mass flow rate and the excess air coefficient on the heating performance of the target soil. The control strategy is explored to accomplish the heating process as expected with constant temperature change rate or constant volumetric water content change rate at different phases by adapting the natural gas flow. The results demonstrate that the heating plan can be realized within the scheduled 36 days and the total natural gas consumption can be reduced by 24% (1487 kg) compared to that of the open-loop reference condition, which may be widely applied for other thermal remediation systems of the polluted soil.


Introduction
Over the past few decades, soil contamination has become an increasingly prominent concern owing to accelerated industrialization and urbanization in China and all over the world [1][2][3]. This phenomenon is particularly severe in the vicinity of heavy industrial areas such as the coking plants, chemical plants and steelworks [4,5]. Among all kinds of hazardous substances released by industrial activities, heavy metals make a significant contribution to soil pollution [6] which are covert, persistent and irreversible [7]. Unlike organic pollutants, which may degrade to less harmful components by means of biological or chemical approaches, heavy metals are not degradable by natural processes. What is worse, this kind of pollution can easily enter into the food chain and also seriously threatens the health and well-being of animals and human beings [8,9]. Due to the severe damage of soil pollution and the severity of the situation, it is urgent that appropriate technologies of contaminated soil treatment are developed.
In situ thermal desorption (ISTD) is one of the remediation technologies for polluted soil and has attracted great research interest in recent years. This method uses fluid injection or heating elements to heat the contaminated site at sufficient temperature without excavation. As a result, most pollutants are separated from soil and collected by the aboveground off-gas treatment units for standardized discharge [10]. Compared with ex situ thermal desorption (ESTD), ISTD has advantages scale apparatus with a small portion of the soil excavated. For field experiments, it is extremely difficult to monitor the moisture transfer process within the heated soil as well as between the heated soil and the surrounding soil along the treatment duration using existing technologies. For laboratory experiments, likewise, owing to that the target recovery soil is excavated and put into a laboratory appliance individually, the moisture transfer between the heated site and the surrounding site during the heating process cannot be simulated veritably and accurately. Due to the disadvantages of the aforementioned experimental research, a lumped-parameter method is employed in the present study in order to investigate the influence of moisture migration on the heating performance of the soil. Furthermore, other correlative studies [31][32][33] by our team using a numerical simulation method or a lumped-parameter method have been conducted and published online.
In this paper, a planned heating control strategy of a natural gas thermal desorption system (GTDS) is proposed aiming to effectively heat the contaminated soil within the specified heating time with less energy consumption. Dynamic heat-transfer modeling for the GTDS on the basis of the lumped-parameter method is established. Simulation analyses are conducted to investigate the impacts of the natural gas mass flow rate and the excess air coefficient on the heating characteristics of the soil. A scheduled heating control strategy based on a group of proportional-integral-derivative (PID) controllers is designed to manipulate the temperature change rate or volumetric water content change rate of the soil by adjusting the natural gas mass flow rate, and the remediation time ultimately. Confronting various disturbances in excess air coefficient and external rainfall, the GTDS under this innovative control strategy can attain target heating temperature in terms of the evidently decreased natural gas consumption and the scheduled operational period with excellent control ability. The results proved that the carefully designed GTDS may have a broad prospect of application for contaminated soil remediation.

General Conception of Gas Thermal Desorption System (GTDS)
The schematic diagram of the GTDS for soil pollution is illustrated in Figure 1. It consists of the heating unit, extraction unit, ground insulation unit and off-gas treatment unit. The heating unit includes a burner fixed on the top of the well, a thermal well inserted vertically into the soil and a combustion blower placed aboveground at the outlet of the well. The natural gas and air are inhaled by the combustion blower into the combustor and burn along with the high temperature flue gas generated. The hot flue gas is injected downward into the inner pipe of the thermal well and then enters the annular space between the inner pipe and outer pipe then, finally, flows upward along the outer pipe and away from the outlet of the outer pipe. During the flow process of flue gas in the thermal well, heat is transferred to the restoration soil which is in contact with the outer pipe of thermal well. Provided there is high temperature for a long period, the volatile contaminants in the soil can vaporize and then be drawn by the extraction blower into the extraction well. This gaseous product is transported to the off-gas treatment unit to be purified for standardized discharge. Additionally, an insulating layer consisting of cellular concrete, insulating brick and rubble with certain thickness individually is laid on the surface of the soil in order to reduce heat loss and prevent pollutant diffusion during operation. Detailed geometric characteristics of the GTDS are provided in Table A1 in Appendix A.
A large number of experimental research works [23,34] have revealed that during the remediation process, the temperature of the heated soil goes through three phases, namely, heat-up phase, boiling phase and superheating phase, which is depicted in Figure 2.
In the heat-up phase, the soil, consisting of solid particles as well as liquid water and vaporous water in the pores, is heated from initial temperature t s0 to water's boiling point t s1 . There exist inappreciable evaporation of liquid water and moisture migration from unheated sites in this stage. In the boiling phase, liquid water evaporates along with a constant temperature of the soil. At this Energies 2020, 13, 642 4 of 28 point, all the heat generated from the thermal well is used to vaporize liquid water. With evaporation of liquid water in the heated soil, the amount of moisture migration from unheated sites increases gradually. When all the liquid water has been vaporized, the target soil enters into the next phase. In the superheating phase, the dry soil can be heated to the final target temperature, meanwhile, the fluid water permeated from unheated zones will evaporate immediately on account of the extremely high temperature. Furthermore, the contaminations can vaporize with the assistant of high temperature for further treatment. A large number of experimental research works [23,34] have revealed that during the remediation process, the temperature of the heated soil goes through three phases, namely, heat-up phase, boiling phase and superheating phase, which is depicted in Figure 2. In the heat-up phase, the soil, consisting of solid particles as well as liquid water and vaporous water in the pores, is heated from initial temperature ts0 to water's boiling point ts1. There exist inappreciable evaporation of liquid water and moisture migration from unheated sites in this stage. In the boiling phase, liquid water evaporates along with a constant temperature of the soil. At this  A large number of experimental research works [23,34] have revealed that during the remediation process, the temperature of the heated soil goes through three phases, namely, heat-up phase, boiling phase and superheating phase, which is depicted in Figure 2. In the heat-up phase, the soil, consisting of solid particles as well as liquid water and vaporous water in the pores, is heated from initial temperature ts0 to water's boiling point ts1. There exist inappreciable evaporation of liquid water and moisture migration from unheated sites in this stage. In the boiling phase, liquid water evaporates along with a constant temperature of the soil. At this point, all the heat generated from the thermal well is used to vaporize liquid water. With evaporation of liquid water in the heated soil, the amount of moisture migration from unheated sites increases gradually. When all the liquid water has been vaporized, the target soil enters into the next phase. In the superheating phase, the dry soil can be heated to the final target temperature, meanwhile, the Note that in the reminder of this paper, the three phases of heat-up phase, boiling phase and superheating phase are defined as Phase-1, Phase-2 and Phase-3 separately for the sake of simplicity. Considering that in the actual construction, the operational period of remediation is generally uncontrollable beforehand but is extremely important not only for the remediation effect but also for the construction progress. To address this problem, a planned heating control strategy based on a group of PID controllers for the GTDS is put forward, as shown in Figure 3. As stated above in Figure 2, the volumetric water content (VWC) in Phase-1 varies little, the temperature during Phase-2 is constant at t s1 , and the VWC at Phase-3 is fixed to 0. Taking this into account, temperature variation rate of Phase-1, VWC change rate of Phase-2 and temperature rise rate of Phase-3 are selected as the target variables, moreover, mass flow rate of natural gas is set as the control variable.

Mathematical Modeling
In this section, the working principle of the GTDS is provided in terms of the mathematics and thermodynamics including the lumped-parameter model and the control model, which may contribute to a better understanding of the function and dynamic characteristics of the system.

Lumped-Parameter Model
As illustrated in Figure 4a, a dynamic heat-transfer model based on the lumped-parameter method for the GTDS is proposed to investigate the thermal performance of the system in detail. In accordance with layout of the GTDS in Figure 1, the burner, the inner pipe of the heat well, the outer pipe of the heat well and the heated soil can be treated as four lumped-parameter nodes respectively, which will be fully described afterwards in this section. Overall, the proposed control strategy comprised a planning layer and a control layer. The fundamental control concept is that by adjusting the natural gas flow, the target soil can be heated with constant temperature rise rate or constant VWC change rate in different stages and attain the objective temperature finally within the scheduled days. To be more specific, the overall heating period and the individual heating time of each phase are firstly determined according to the construction plan. According to both the heating time of each stage and the variation of each objective variable, the theoretical value is further calculated. The calculated value and actual value of each objective variable at different stages are respectively regarded as two inputs of the corresponding PID controllers. Based on the errors of the objective variables, the natural gas can be manipulated with corresponding mass flow rate and transported to the burner as the energy source of the GTDS.

Mathematical Modeling
In this section, the working principle of the GTDS is provided in terms of the mathematics and thermodynamics including the lumped-parameter model and the control model, which may contribute to a better understanding of the function and dynamic characteristics of the system.

Lumped-Parameter Model
As illustrated in Figure 4a, a dynamic heat-transfer model based on the lumped-parameter method for the GTDS is proposed to investigate the thermal performance of the system in detail. In accordance with layout of the GTDS in Figure 1, the burner, the inner pipe of the heat well, the outer pipe of the heat well and the heated soil can be treated as four lumped-parameter nodes respectively, which will be fully described afterwards in this section.

Model of Burner
The schematic diagram of the burner is depicted in Figure 4b. Based on the energy conservation principle, the temperature dynamics of the burner can be expressed by [35,36]: where mb is the mass of the burner; tb and t0 are the temperature of the burner and the ambient temperature; Gg, Ga and Gf are the mass flow rates of the natural gas, the air and the flue gas respectively; Lg is the volume flow rate of the natural gas; QL is low heating value of the natural gas; Rb is thermal resistance between the burner wall and burner crick; hg, ha and hf are the enthalpies of the natural gas, the air and the flue gas which can be written as: where tb, tg and ta are the temperatures of the burner, the natural gas and the air respectively; cg, ca and cf are the specific heat of the natural gas, the air and the high temperature flue gas respectively. Thereby, cf can be calculated by:

Model of Burner
The schematic diagram of the burner is depicted in Figure 4b. Based on the energy conservation principle, the temperature dynamics of the burner can be expressed by [35,36]: where m b is the mass of the burner; t b and t 0 are the temperature of the burner and the ambient temperature; G g , G a and G f are the mass flow rates of the natural gas, the air and the flue gas respectively; L g is the volume flow rate of the natural gas; Q L is low heating value of the natural gas; R b is thermal Energies 2020, 13, 642 7 of 28 resistance between the burner wall and burner crick; h g , h a and h f are the enthalpies of the natural gas, the air and the flue gas which can be written as: where t b , t g and t a are the temperatures of the burner, the natural gas and the air respectively; c g , c a and c f are the specific heat of the natural gas, the air and the high temperature flue gas respectively. Thereby, c f can be calculated by: where g i is mass fraction of the flue gas, c fi is mass specific heat of components of the flue gas; the output flue gas temperature from the burner t f can be expressed by: where η b is the heat exchange efficiency of the burner.
In addition, G a and G f are given by Equations (5) and (6) respectively, where V 0 is theoretical volume of the air for complete combustion of the natural gas, α is excess air coefficient.
2. Model of Thermal Well Figure 4c illustrates individual temperature points of the heat well. The temperature dynamics of the inner pipe and outer pipe of the thermal well, which are based on the energy conservation principle, are represented by Equations (7)-(8) [31,37]: (m W2 c W2 ) where m W1 and m W2 are the mass of the inner pipe and outer pipe; t W1 , t W2 and t s are the temperatures in Celsius of the inner pipe, the outer pipe and the heated soil respectively; T W1 and T W2 are Kelvin temperatures of the inner pipe and outer pipe; t fi1 and t fi2 are the input temperatures of the flue gas that flows into the inner pipe and the outer pipe; c W1 and c W2 are the specific heat of the inner pipe and outer pipe; η W is the heat-exchange efficiency of the thermal well; R Ws is thermal resistance between the outer pipe and the soil; C 0 is the black-body radiation coefficient; A 1 and A 2 are surface area of the inner pipe and outer pipe; ε 1 and ε 2 are the emissivity of the inner pipe and outer pipe; ε a is the system emissivity which is defined in Equation (9) [38].
Equations (10) and (11) present the output flue gas temperatures from the inner pipe (t fo1 ) and that from the outer pipe (t fo2 ) respectively. As clearly depicted in Figure 4c, t fo1 is actually the input flue gas temperature of the outer pipe (t fi2 ), which is expressed by Equation (12).
Energies 2020, 13, 642 8 of 28 3. Model of Heated Soil For simplicity, some assumptions are made as follows: (1) The soil is homogeneous porous medium and there is no chemical interaction in the soil.
(2) Solid, liquid, and gas phases are separately continuous in unsaturated soil.
(3) Natural convection of fluids in the soil satisfies Darcy's law. (4) The target site is considered as one of the unit blocks and each block is heated by a thermal well individually. There is no heat exchange between the heated site and the surrounding sites at the boundary. (5) There is no moisture transfer between the target heated site and the surrounding heated sites, moreover, the moisture transfer mainly occurs between the heated soil and the underlying unheated soil. (6) The influence of gravity on the soil moisture transfer in the unsaturated soil is assumed to be negligible. (7) The pressure inside the soil is evenly distributed. (8) The gases in the soil are assumed to be ideal.
(a) Liquid flow model Based on the mass conservation principle, the VWC variation Equation [39] can be written as: where θ s is the VWC of the soil, ρ w is density of liquid water, V s is volume of the soil, m wv represents mass of vaporized liquid water, and m wi denotes mass of liquid water flowing vertically into the heated soil [40] which is given by: where A cs is the migration area of liquid water in the vertical direction, J wi is the liquid water migration flux density. A lot of research [41][42][43][44] has shown that water seepage in the unsaturated soil is primarily caused by the temperature gradient and the moisture gradient. The liquid water migration flux density J wi can be expressed by Equation (15) in line with Assumptions (5)-(6), where t s is temperature of the heated soil, D wt and D wθ are the mass diffusivities of the liquid water in soil caused by temperature gradient and moisture gradient respectively; l is the transfer distance of the liquid water in the vertical direction, K w is the hydraulic conductivity of the soil [45,46] which is defined by: ε is the void ratio of the soil, K ws is saturated hydraulic conductivity, b is the characteristic coefficient of the soil.
(b) Vapor flow model Likewise, on the basis of mass conservation principle, the mass equation of vapor can be written as: Energies 2020, 13, 642 where V vs is vapor volume in the soil, ρ v is density of vapor water, m vo is mass of vapor water extracted into the extraction well which is defined by: where N e is the power of the extraction blower, P eo is the outlet pressure of the extraction blower, R v is the gas constant of vapor water, T v is Kelvin temperature of the vapor water. Based on Assumption (7), provided that the volume of vapor in the soil pores V vs and that in the extraction well V v_e are considered as a whole, the pressure variation equation can be derived from Equations (17) and (18) that: (c) Heat-Flow Model As presented in Figure 5, the temperature dynamics of the soil [31,33] can be expressed by: where Q W is heating power of the thermal well; Q up and Q down are heat leakage between the heated site and the insulating layer, as well as that between the heated site and the underlying unheated soil; Q W , Q up and Q down can be calculated by: Energies 2020, 13, 642 10 of 29 where QW is heating power of the thermal well; Qup and Qdown are heat leakage between the heated site and the insulating layer, as well as that between the heated site and the underlying unheated soil; QW, Qup and Qdown can be calculated by: Qw and Qv represent the heat flux of liquid water migration and vapor water extraction respectively which can be given by: where hw and hv are specific enthalpy of liquid water and vaporous water, mwi and mvo can be calculated by Equations (14) and (18). Qeva indicates the energy absorbed by liquid water evaporation which is defined by: where H is latent heat of evaporation.

Fluid-Flow Model
As shown in Figure 6, the overall flow resistance network of the gaseous fluids for the GTDS is a non-linear parallel and series flow network which can be expressed in detail by the following equations. Specifically, the fluid flow network comprises of the inlet natural gas pipe and inlet air pipe which are connected in parallel, as well as the burner wall, the heat well inner pipe, the heat well outer pipe and the flue gas outlet pipe which are associated in series. When the gaseous fluids such as natural gas, air and flue gas flow through the above channels, there exist flow resistances Q w and Q v represent the heat flux of liquid water migration and vapor water extraction respectively which can be given by: where h w and h v are specific enthalpy of liquid water and vaporous water, m wi and m vo can be calculated by Equations (14) and (18).
Q eva indicates the energy absorbed by liquid water evaporation which is defined by: where H is latent heat of evaporation.

Fluid-Flow Model
As shown in Figure 6, the overall flow resistance network of the gaseous fluids for the GTDS is a non-linear parallel and series flow network which can be expressed in detail by the following equations. Specifically, the fluid flow network comprises of the inlet natural gas pipe and inlet air pipe which are connected in parallel, as well as the burner wall, the heat well inner pipe, the heat well outer pipe and the flue gas outlet pipe which are associated in series. When the gaseous fluids such as natural gas, air and flue gas flow through the above channels, there exist flow resistances correspondingly which can be expressed by Equation (24) [47], where µ is the fluid viscosity, l is the channel length, ξ is the equivalent coefficient of channel length considering the bend and corner [48].
Energies 2020, 13, 642 11 of 29 Given the flow resistance of the inlet natural gas pipe Rg and that of the inlet air pipe Ra calculated from Equation (24), the parallel flow resistance Rpar [49,50] can be presented by Equation (25), where ϕg and ϕa characterize the valve opening of the inlet natural gas pipe and air pipe. After that, the overall flow resistance Rnet is proposed in Equation (26) where Rb, Rwi, Rwo and Rpo are the flow resistances corresponding to each series branches as expressed in Figure 6. .
For a given pressure difference ∆Pc defined in Equation (27), where Pci and Pco denote the inlet pressure and outlet pressure of the combustion blower respectively, the mass flow rate of flue gas through the fluid loop Gf can be written as Equation (28).
According to Equations (5) and (6), the natural gas mass flow rate Gg and the air mass flow rate Ga can be separately expressed by Equations (29) and (30), correspondingly. Furthermore, the natural gas inlet pressure Pg and the air inlet pressure Pa can be obtained by Equations (31) and (32), where Pb is inlet pressure of the burner. Given the flow resistance of the inlet natural gas pipe R g and that of the inlet air pipe R a calculated from Equation (24), the parallel flow resistance R par [49,50] can be presented by Equation (25), where φ g and φ a characterize the valve opening of the inlet natural gas pipe and air pipe. After that, the overall flow resistance R net is proposed in Equation (26) where R b , R wi , R wo and R po are the flow resistances corresponding to each series branches as expressed in Figure 6.
For a given pressure difference ∆P c defined in Equation (27), where P ci and P co denote the inlet pressure and outlet pressure of the combustion blower respectively, the mass flow rate of flue gas through the fluid loop G f can be written as Equation (28).
According to Equations (5) and (6), the natural gas mass flow rate G g and the air mass flow rate G a can be separately expressed by Equations (29) and (30), correspondingly. Furthermore, the natural gas inlet pressure P g and the air inlet pressure P a can be obtained by Equations (31) and (32), where P b is inlet pressure of the burner.

Planned Heating Control Model
The block diagram of the GTDS with planned heating control is shown in Figure 7. On the basis of conventional PID control, the amount of the natural gas is adjusted according to the difference of temperature variation rate or VWC variation rate at different stages, and finally achieve the heating plan in a specific number of days with less energy consumption. Whatever phase the soil is heated to, the temperature variation and the VWC variation can be calculated by: where n is the step size.
Energies 2020, 13, 642 12 of 29 where n is the step size.
As inputs of the three PID controllers, the control error e1, e2 and e3 are described by: where ' 1 ts  With the assistance of the merge module, the mass flow rate of the natural gas (Gg) which is the energy input of the burner can be given by Equation (35), where Gg1, Gg2, Gg3 are the mass flow rate of the natural gas at each stage respectively.

Simulation Condition Settings
Related parameter determinations of the GTDS and the original state of the system are summarized in Table A2 in Appendix B. The mass flow rate of the natural gas (Gg) and excess air coefficient (α) are two primary input variables which will be changed and controlled throughout the simulation process. As the base line for transient performance and control effect analyses, the initial Gg is 0.989 × 10 −3 kg/s, and the initial α is 2.2.
Matlab R2017a was used as the simulation software. The flowchart of the simulation process is illustrated in Figure 8 which can be described in detail as follows: (1) the program is initialized at the beginning in accordance with the parameters tabulated in Table A2. (2)   As inputs of the three PID controllers, the control error e 1 , e 2 and e 3 are described by: where ν ts1 , ν θs2 and ν ts3 are desired temperature change rate in Phase-1, the wanted VWC change rate in Phase-2 and the desired temperature change rate in Phase-3 respectively; v ts1 , v ts2 and v ts3 are actual values correspondingly.
With the assistance of the merge module, the mass flow rate of the natural gas (G g ) which is the energy input of the burner can be given by Equation (35), where G g1 , G g2 , G g3 are the mass flow rate of the natural gas at each stage respectively.

Simulation Condition Settings
Related parameter determinations of the GTDS and the original state of the system are summarized in Table A2 in Appendix B. The mass flow rate of the natural gas (G g ) and excess air coefficient (α) are two primary input variables which will be changed and controlled throughout the simulation process. As the base line for transient performance and control effect analyses, the initial G g is 0.989 × 10 −3 kg/s, and the initial α is 2.2.
Matlab R2017a was used as the simulation software. The flowchart of the simulation process is illustrated in Figure 8 which can be described in detail as follows: (1) the program is initialized at the beginning in accordance with the parameters tabulated in Table A2. (2) The simulation time and the calculation step are set. (3) Initial parameters are provided. (4) The desired values of ν ts1 , ν θs2 and ν ts3 are derived. (5) The actual values of v ts1 , v ts2 and v ts3 can be obtained by Equations (1)-(23) and Equation (33). (6) During the closed-loop control, the control errors e 1 , e 2 and e 3 can be calculated by Equation (34) firstly; and then G g can be adjusted according to Equation (35). (7) After all the above procedures, two judgments need to be made. The first judgment is whether the control objective is convergence. If no, the simulation process will be updated by the new control parameters. If the first judgment is yes, the simulation continues to the second judgment, which is whether the simulation will continue. If the second judgment is yes, the simulation will jump to step (2). If the second judgment is no, the simulation will end.
Energies 2020, 13, 642 13 of 29 Equation (33). (6) During the closed-loop control, the control errors e1, e2 and e3 can be calculated by Equation (34) firstly; and then Gg can be adjusted according to Equation (35). (7) After all the above procedures, two judgments need to be made. The first judgment is whether the control objective is convergence. If no, the simulation process will be updated by the new control parameters. If the first judgment is yes, the simulation continues to the second judgment, which is whether the simulation will continue. If the second judgment is yes, the simulation will jump to step (2). If the second judgment is no, the simulation will end. For investigating the effects of various input variables on the heating performance of the target soil, two simulation conditions are organized and listed in Table 1. Two simulation conditions are organized and listed in Table 1. To be specific, various step disturbances in the input Gg and α take place to demonstrate the influence of the natural gas mass flow rate and the excess air coefficient upon the soil's heating characteristics.   For investigating the effects of various input variables on the heating performance of the target soil, two simulation conditions are organized and listed in Table 1. Two simulation conditions are  Table 1. To be specific, various step disturbances in the input G g and α take place to demonstrate the influence of the natural gas mass flow rate and the excess air coefficient upon the soil's heating characteristics. Given a fixed temperature variation or VWC variation, the desired temperature change rate (ν ts ) or desired VWC change rate (ν θs ) vary with the heating time at different phases. For the purpose of investigating the influence of ν ts or ν θs on the system performances, simulation conditions for each phase are offered in Table 2. It is worth noting that the initial temperature and the objective temperature of soil are assumed to be 20 • C and 525 • C, respectively, which leads to the temperature change of Phase-1 being 80 • C and that of Phase-3 being 425 • C. Moreover, considering an original VWC of 0.25 m 3 /m 3 and the VWC increment of 0.0005 m 3 /m 3 in Phase-1, the VWC variation of Phase-2 is −0.2505 m 3 /m 3 . Determinations of relevant PID parameters are summarized in Table 3.   Table 3. Parameters of the simulated proportional-integral-derivative (PID) controllers.

Controllers Parameters
In addition, to investigate the performance in rejecting disturbances of the controller, two cases are arranged and listed in Table 4. Step disturbance G wi 100% step increase lasting for 4 days in each of the three phases

Open-Loop Dynamic Characteristics
The focus of this section is to investigate the effects of the natural gas mass flow rate G g and excess air coefficient α on the temperature characteristic and moisture characteristic of the heated soil. Due to that the soil heating process is a slowly varying process, the sampling step is set as 3600 s (1 h) for simulation. Figure 9 presents the temperature and VWC characteristics of the target soil under various step disturbances in G g corresponding to the first simulation condition in Table 1. Note that the G g of 0.989 × 10 −3 kg/s is considered as a baseline working condition with other four conditions. In Figure 9a, take the baseline condition (0.989 × 10 −3 kg/s) which is represented by the blue for instance, t s rises from 20 • C to 100 • C within the first 8.5 days (Phase-1), and then is fixed at 100 • C in the next 24.8 days (Phase-2), in Phase-3, it goes up rapidly from 100 • C and settles at 532 • C ultimately in the last 39 days. For other cases, t s presents the same trend as well. The variation trend of the soil temperature is highly consistent with the research results in related literatures [23,34].

Disturbance of Natural Gas Mass Flow Rate
The variation of t s can be further elaborated by v ts which is shown in Figure 9b. For the baseline case in Phase-1, v ts decreases slowly within the range from 0.42 • C /h to 0.36 • C /h; in Phase-2, v ts keeps 0 all the time; in Phase-3, v ts raises sharply to the maximum of 2.03 • C /h firstly, then declines gradually and attains stabilization of nearly 0 • C/h ultimately. For other conditions, the overall trend is the same as the baseline condition. Figure 9c provides the VWC (θ s ) variation under various G g . When G g is 0.989 × 10 −3 kg/s, θ s exhibits a slight increase from 0.25 m 3 /m 3 to 0.2505 m 3 /m 3 at Phase-1; in Phase-2, θ s declines considerably from 0.2505 m 3 /m 3 to 0 m 3 /m 3 owing to the large evaporation amount of liquid water; during Phase-3, θ s remains unchanged of 0 m 3 /m 3 . For further explaining, Figure 9d reflects the variation of v θs under different G g . In accordance with θ w , v θs is nearly 0 in Phase-1 and remains 0 during Phase-3, but presents a slow rise in Phase-2. during Phase-3, but presents a slow rise in Phase-2.
The individual heating time of three phases τ1, τ2 and τ3 as well as the total heating period τs at different Gg are plotted in Figure 9e. With the increase in Gg, τ1 decreases from 17.6 days to 6.25 days, τ2 descends from 55.2 days to 17.3 days, τ3 drops from 47.9 days to 32.8 days. In addition, the sum of τ1, τ2 and τ3, τs keeps the downward tendency as well. The stable temperature of the soil (ts_sta) at various Gg is depicted in Figure 9f. It can be viewed that ts_sta raises exponentially from 358 °C to 619 °C with the increase in Gg from −50% step-disturbance to +50% step-disturbance. The above observations demonstrate that the natural gas consumption has a decisive effect on the heating performance of the GTDS absolutely. As for the baseline condition (Gg = 0.989 × 10 −3 kg/s), the heating period is 72.3 days and the natural gas consumption is 6182 kg. More natural gas applied to the GTDS leads to a higher stable temperature and a shorter remediation time of the target site. However, the higher the Gg, the more energy the system will consume, which is not preferable. Aiming to balance the energy consumption and the heating duration which are two dominating parameters in this paper, control strategy analyses in different stages will be conducted in Section 4.1.

Disturbance of Excess Air Coefficient
The characteristics of soil under different disturbances in α in line with the second simulation condition in Table 1 are expressed in Figure 10. Figure 10a,b present the variations of ts and vts with various α. Taking the case of 1.87 for example, ts increases with a slightly decreased vts from in Phase-1; and then both ts and v are unchanged all the while in Phase-2; finally, ts rises rapidly with a sharply The individual heating time of three phases τ 1 , τ 2 and τ 3 as well as the total heating period τ s at different G g are plotted in Figure 9e. With the increase in G g , τ 1 decreases from 17.6 days to 6.25 days, τ 2 descends from 55.2 days to 17.3 days, τ 3 drops from 47.9 days to 32.8 days. In addition, the sum of τ 1 , τ 2 and τ 3 , τ s keeps the downward tendency as well. The stable temperature of the soil (t s_sta ) at various G g is depicted in Figure 9f. It can be viewed that t s_sta raises exponentially from 358 • C to 619 • C with the increase in G g from −50% step-disturbance to +50% step-disturbance. The above observations demonstrate that the natural gas consumption has a decisive effect on the heating performance of the GTDS absolutely. As for the baseline condition (G g = 0.989 × 10 −3 kg/s), the heating period is 72.3 days and the natural gas consumption is 6182 kg. More natural gas applied to the GTDS leads to a higher stable temperature and a shorter remediation time of the target site. However, the higher the G g , the more energy the system will consume, which is not preferable. Aiming to balance the energy consumption and the heating duration which are two dominating parameters in this paper, control strategy analyses in different stages will be conducted in Section 4.1.

Disturbance of Excess Air Coefficient
The characteristics of soil under different disturbances in α in line with the second simulation condition in Table 1 are expressed in Figure 10. Figure 10a,b present the variations of t s and v ts with various α. Taking the case of 1.87 for example, t s increases with a slightly decreased v ts from in Phase-1; and then both t s and v ts are unchanged all the while in Phase-2; finally, t s rises rapidly with a sharply declined v ts in Phase-3. The duration of three phases are 8.3 days, 23.9 days and 40.96 days, respectively. The VWC (θs) variation and VWC change rate (vθs) variation under different α are graphed in Figure 10c,d respectively. For all the five conditions, the overall variant trend of θs and vθs are similar to that in Figure 9c,d in Section 3.1. Note that vθs is nearly close to 0 and θs exhibits a slight increase in Phase-1; θs and vθs remains 0 in Phase-3; whereas in Phase-2, θs goes down to 0 m 3 /m 3 gradually and the corresponding vθs shows a slow rise. Furthermore, for the second phase, the greater the α is, the smaller the absolute value of vθs, and the longer it takes for θs to decline to 0 m 3 /m 3 .
As revealed in Figure 10e, when confronting disturbances in α ranging from 1.54 to 2.86, τ1 presents a slight rise from 8.2 days to 9.3 days, τ2 exhibits an obvious increase from 23.3 days to 27 The VWC (θ s ) variation and VWC change rate (v θs ) variation under different α are graphed in Figure 10c,d respectively. For all the five conditions, the overall variant trend of θ s and v θs are similar to that in Figure 9c,d in Section 3.1. Note that v θs is nearly close to 0 and θ s exhibits a slight increase in Phase-1; θ s and v θs remains 0 in Phase-3; whereas in Phase-2, θ s goes down to 0 m 3 /m 3 gradually and the corresponding v θs shows a slow rise. Furthermore, for the second phase, the greater the α is, the smaller the absolute value of v θs , and the longer it takes for θ s to decline to 0 m 3 /m 3 .
As revealed in Figure 10e, when confronting disturbances in α ranging from 1.54 to 2.86, τ 1 presents a slight rise from 8.2 days to 9.3 days, τ 2 exhibits an obvious increase from 23.3 days to 27 days, whereas, τ 3 decreases considerably from 44.8 days to 30.8 days. The combination of τ 1 , τ 1 and τ 3 yields that τ s declines by 9.3 days accordingly. Figure 10f reflects that t s_sta decreases linearly at a rate of −125 approximately with the increase in α.
The analyses stated above suggest that the excess air coefficient exerts a significant influence upon the soil's temperature and VWC characteristics. Provided that the mass flow rate of natural gas is constant and the excess air coefficient is greater than 1.05, less air brought to the burner results in a higher stable temperature and a shorter heating period. This can be explained by the fact that proper increase of excess air coefficient can reduce the heat loss of chemical incomplete combustion; too high an excess air coefficient, however, will reduce the furnace temperature and yield incomplete combustion.

Closed-Loop Control Effect Analyses
In this section, the GTDS with the well-designed planned heating controller is firstly discussed in terms of various temperature change rates or VWC change rates in different phases, and then responses to disturbances in the excess air coefficient and liquid water migration quantity for closed-loop simulation, which is based on a group of representative theoretical values of objective variables (ν ts1 , ν θs2 and ν ts3 ).

Closed-Loop Control Performance of Phase-1
The closed-loop control performances of the GTDS under different ν ts1 are expressed in Figure 11 corresponding to the first simulation condition in Table 2. It can be seen from Figure 11a that when confronting the increasing τ 1 from 4 days to 10 days, t s raises linearly at a constant rate under the control of PID-1. As mentioned in Section 3.1, v ts1 will decline without closed-loop control. As shown in Figure 11b, for all the temperature change rate variations, v ts1 achieve the objectives under closed-loop control within the separate planned time. However, given the same control parameters of PID-1 displayed in Table 3, there are several specific differences in the control effects among these v ts1 variation conditions. Generally speaking, greater theoretical value of the temperature rise rate results in smaller overshoot, shorter settling time and larger steady-state error. Take the ν ts1 of 0.5556 • C/h for instance, the overshoot, settling time, and steady-state error are calculated to be 25.36% (0.14 • C/h), 15 h and 0.28% (0.0016 • C/h), respectively.
As depicted in Figure 11c, G g presents a gradual increasing trend. Furthermore, the greater the ν ts1 is, the more obviously the G g increases. According to the above simulation results, the heating duration (τ 1 ) and the total consumption of the natural gas (M g1 ) of Phase-1 under different temperature rise rates is plotted in Figure 11d. It can be seen that τ 1 declines gradually with the increasing ν ts1 and the relationship between them is approximately inversely proportional. In addition, M g1 increases slowly when ν ts1 is less than 0.556 • C/h and M g1 increases rapidly when ν ts1 is greater than 0.556 • C/h. In summary, the heating plan in the first stage can be realized within specified days by adjusting the change rate of soil's temperature. Increasing the temperature rise rate can significantly shorten the heating time of the first stage, whereas, at the cost of increasing the consumption of the natural gas. Aiming to balance the construction duration and the energy consumption, 0.556 °C/h is set as the expected temperature change rate of Phase-1 in the following closed-loop simulation of Phase-2 and Phase-3.

Dynamic Control Performance of Phase-2
On the premise that the theoretical temperature variation rate of Phase-1 is 0.556 °C/h, PID-2 takes effect with various VWC change rates and the simulation results are presented in Figure 12. As stated in Figure 9c in Section 3.1, at the end of Phase-1, θs goes up to 0.2505 m 3 /m 3 which is also the starting value of the VWC of Phase-2. As shown in Figure 12a, θs descends linearly from 0.2505 m 3 /m 3 to 0 m 3 /m 3 at a unchanged rate (vθs2) with the increasing heating time of Phase-2 (τ2) from 10 days to 30 days with the increment of 5 days.
It can be seen from Figure 12b that, for all the VWC change rate variations of Phase-2, vθs2 settles to the target value within the respective desired time under closed-loop control. To be specific, under the same control parameters of PID-2 tabulated in Table 3, the greater the absolute value of vθs2 is, the greater the overshoot, the shorter the settling time, and the larger the steady-state error that will be obtained. When ' 2 s θ ν is −0.697 × 10 −3 m 3 /(m 3 ·h), the overshoot and settling time are, respectively, 9.32% (0.65 × 10 −4 m 3 /(m 3 ·h) and 9 h, and the steady-state error is infinitesimally small and can be ignored. In summary, the heating plan in the first stage can be realized within specified days by adjusting the change rate of soil's temperature. Increasing the temperature rise rate can significantly shorten the heating time of the first stage, whereas, at the cost of increasing the consumption of the natural gas. Aiming to balance the construction duration and the energy consumption, 0.556 • C/h is set as the expected temperature change rate of Phase-1 in the following closed-loop simulation of Phase-2 and Phase-3.

Dynamic Control Performance of Phase-2
On the premise that the theoretical temperature variation rate of Phase-1 is 0.556 • C/h, PID-2 takes effect with various VWC change rates and the simulation results are presented in Figure 12. As stated in Figure 9c in Section 3.1, at the end of Phase-1, θ s goes up to 0.2505 m 3 /m 3 which is also the starting value of the VWC of Phase-2. As shown in Figure 12a, θ s descends linearly from 0.2505 m 3 /m 3 to 0 m 3 /m 3 at a unchanged rate (v θs2 ) with the increasing heating time of Phase-2 (τ 2 ) from 10 days to 30 days with the increment of 5 days.
It can be seen from Figure 12b that, for all the VWC change rate variations of Phase-2, v θs2 settles to the target value within the respective desired time under closed-loop control. To be specific, under the same control parameters of PID-2 tabulated in Table 3, the greater the absolute value of v θs2 is, the greater the overshoot, the shorter the settling time, and the larger the steady-state error that will be obtained. When ν θs2 is −0.697 × 10 −3 m 3 /(m 3 ·h), the overshoot and settling time are, respectively, 9.32% (0.65 × 10 −4 m 3 /(m 3 ·h) and 9 h, and the steady-state error is infinitesimally small and can be ignored. It can be concluded from the above results that the boiling phase of the target soil can be accomplished as planned via the control of the VWC change rate. Higher absolute value of the VWC variation rate of the second stage leads to shorter heating time, which is beneficial, as well as more consumption of the natural gas, which is unfavorable. Therefore, −0.697 × 10 −3 m 3 /(m 3 ·h) is chosen as the desired VWC change rate of Phase-2 for the closed-loop simulation of Phase-3.

Characteristics of Phase-3 under Closed-Loop Control
Provided a fixed temperature rise rate of 0.556 ℃/h for Phase-1 and an invariable VWC change rate of −0.697 × 10 −3 m 3 /(m 3 ·h) for Phase-2, the closed-loop simulation of the GTDS under the control of PID-3 is conducted in line with the third simulation condition in Table 2, and the relevant results are illustrated in Figure 13.
It can be observed from Figure 13a that under the control of PID-3, ts rises linearly from 100 °C to the target 525 °C at a constant rate (vts3) which is offered in Figure 13b. As the objective variable of Confronting ν θs2 varying conditions, the G g in Figure 12c remains almost constant all through Phase-2 after a transient oscillation. Take the ν θs2 of −0.697 × 10 −3 m 3 /(m 3 ·h) for example, after an oscillation of 9 h, the G g settles to be 0.001814 kg/s approximately. In addition, the construction time (τ 2 ) and the mass of the natural gas consumed (M g2 ) are shown graphically in Figure 12d. Considering that the v θw2 is negative, Figure 12d takes |ν θs2 | as the abscissa. With the |ν θs2 | raising, τ 2 exhibits an decrease as expected. Meanwhile, the M g2 increases slowly when |ν θs2 | is less than 0.697 × 10 −3 m 3 /(m 3 ·h) and rises sharply when |ν θs2 | is greater than 0.697 × 10 −3 m 3 /(m 3 ·h).
It can be concluded from the above results that the boiling phase of the target soil can be accomplished as planned via the control of the VWC change rate. Higher absolute value of the VWC variation rate of the second stage leads to shorter heating time, which is beneficial, as well as more consumption of the natural gas, which is unfavorable. Therefore, −0.697 × 10 −3 m 3 /(m 3 ·h) is chosen as the desired VWC change rate of Phase-2 for the closed-loop simulation of Phase-3.

Characteristics of Phase-3 under Closed-Loop Control
Provided a fixed temperature rise rate of 0.556 • C/h for Phase-1 and an invariable VWC change rate of −0.697 × 10 −3 m 3 /(m 3 ·h) for Phase-2, the closed-loop simulation of the GTDS under the control of PID-3 is conducted in line with the third simulation condition in Table 2, and the relevant results are illustrated in Figure 13. steady-state error of vts3 increases slightly and gradually due to the sharp decline of the temperature rise rate during Phase-3 without closed-loop control. As the control variable of PID-3, meanwhile, Gg in Figure 13c presents considerable increase corresponding to all the ' 3 ts ν variation conditions. In addition, the greater the ' 3 ts ν is, the more rapidly the Gg rises. For the maximal ' The heating time (τ3) and the total natural gas consumption (Mg3) of Phase-3 versus temperature rise rate are described in Figure 13d. With the increasing vts3 from 0.59 °C/h to 1.77 °C/h, τ3 varies from 30 days to 10 days accordingly where the relationship between τ3 and vts3 suggests an approximately negative proportional relationship. In the meantime, Mg3 declines in a similar way from 2112 kg to 1470 kg with the climbing vts3.
The observations suggest that increasing the temperature change rate of the third stage can not only shorten the heating time greatly, but also reduce the total natural gas consumption considerably, which is of great significance to the soil remediation of the GTDS.
It is worth noting that according to the analyses stated above, the theoretical values of objective variables at different stages are assumed as follows: (1)   It can be observed from Figure 13a that under the control of PID-3, t s rises linearly from 100 • C to the target 525 • C at a constant rate (v ts3 ) which is offered in Figure 13b. As the objective variable of PID-3, v ts3 achieves the objective with a tiny steady-state error for all the cases. In particular, the steady-state error of v ts3 increases slightly and gradually due to the sharp decline of the temperature rise rate during Phase-3 without closed-loop control. As the control variable of PID-3, meanwhile, G g in Figure 13c presents considerable increase corresponding to all the ν ts3 variation conditions. In addition, the greater the ν ts3 is, the more rapidly the G g rises. For the maximal ν ts3 (1.77 • C/h) condition, G g goes up from 0.00094 kg/s to 0.0034 kg/s within the Phase-3 duration of 10 days after a fluctuation of 12 h.
The heating time (τ 3 ) and the total natural gas consumption (M g3 ) of Phase-3 versus temperature rise rate are described in Figure 13d. With the increasing v ts3 from 0.59 • C/h to 1.77 • C/h, τ 3 varies from 30 days to 10 days accordingly where the relationship between τ 3 and v ts3 suggests an approximately negative proportional relationship. In the meantime, M g3 declines in a similar way from 2112 kg to 1470 kg with the climbing v ts3 .
The observations suggest that increasing the temperature change rate of the third stage can not only shorten the heating time greatly, but also reduce the total natural gas consumption considerably, which is of great significance to the soil remediation of the GTDS.
It is worth noting that according to the analyses stated above, the theoretical values of objective variables at different stages are assumed as follows: (1) the desired temperature change rate of Phase-1 is 0.556 • C/h; (2) the expected VWC variation rate of Phase-2 is −0.697 × 10 −3 m 3 /(m 3 ·h); (3) the target temperature rise rate of Phase-3 is 1.18 • C/h. The overall heating time and the total natural gas consumption thus obtained are 36 days and 4708 kg, respectively. Compared with the open-loop reference condition, the heating time is shortened by 50.32% (36.46 days) and the natural gas consumption is reduced by 24% (1487 kg). On this basis, two kinds of disturbances are brought to the GTDS to examine the closed-loop control effects in rejecting disturbances of the controller.

Excess Air Coefficient Disturbance Response
The system dynamic characteristics response to a periodic disturbance in the input excess air coefficient α to the GTDS are plotted in Figure 14. As revealed in Figure 14a, the periodic disturbance in α is simulated in the form of a square wave pulse with high level of 2.5 and low level of 2.2, respectively. The results show that v ts1 in Figure 14b fluctuates within 4.3% of its steady-state value (the absolute overshoot is 0.0238 • C/h) with the settling time of 9 h, while v θs2 in Figure 14c reaches an equilibrium after 10 s with a 5.8% fluctuation of its steady-state value (the absolute overshoot is 0.4 × 10 −4 m 3 /(m 3 ·h)). Need of special note is that v ts3 in Figure 14b presents a gradually increasing fluctuation along the time. This may be ascribed to the slight decline of the temperature change rate at Phase-3 under the control of PID-3, as stated in Section 4.1.3. Additionally, Figure 14d reflects the fact that G g fluctuates with a slow increasing trend at Phase-1 and attains equilibrium with a 6% fluctuation of 0.00192 kg/s during Phase-2, followed by growing undulation with a rapidly increased central value in Phase-3.
Energies 2020, 13,642 22 of 29 consumption is reduced by 24% (1487 kg). On this basis, two kinds of disturbances are brought to the GTDS to examine the closed-loop control effects in rejecting disturbances of the controller.

Excess Air Coefficient Disturbance Response
The system dynamic characteristics response to a periodic disturbance in the input excess air coefficient α to the GTDS are plotted in Figure 14. As revealed in Figure 14a, the periodic disturbance in α is simulated in the form of a square wave pulse with high level of 2.5 and low level of 2.2, respectively. The results show that vts1 in Figure 14b fluctuates within 4.3% of its steady-state value (the absolute overshoot is 0.0238 °C/h) with the settling time of 9 h, while vθs2 in Figure 14c reaches an equilibrium after 10 s with a 5.8% fluctuation of its steady-state value (the absolute overshoot is 0.4 × 10 −4 m 3 /(m 3 ·h)). Need of special note is that vts3 in Figure 14b presents a gradually increasing fluctuation along the time. This may be ascribed to the slight decline of the temperature change rate at Phase-3 under the control of PID-3, as stated in Section 4.1.3. Additionally, Figure 14d reflects the fact that Gg fluctuates with a slow increasing trend at Phase-1 and attains equilibrium with a 6% fluctuation of 0.00192 kg/s during Phase-2, followed by growing undulation with a rapidly increased central value in Phase-3. Therefore, it can be concluded that when experiencing a periodic disturbance in the excess air coefficient, the well-designed control strategy manifests fast response ability and strong selfadaptability as previously expected. Moreover, the periodic disturbance in the excess air coefficient exerts little influence upon the overall consumption of the natural gas. Therefore, it can be concluded that when experiencing a periodic disturbance in the excess air coefficient, the well-designed control strategy manifests fast response ability and strong self-adaptability as previously expected. Moreover, the periodic disturbance in the excess air coefficient exerts little influence upon the overall consumption of the natural gas.

External Rainfall Disturbance Response
The transient responses of the GTDS to an unfavorable external rainfall disturbance are depicted in Figure 15. From the point of view that the occurrence of rainfall during the remediation will change the amount of liquid water diffusing into the heated soil firstly and affect the heating period and the energy consumption further, G wi with step disturbances is used to simulate the rainfall events. Assuming that the rainfall occurs at each stage and lasts for 4 days, as well as that during periods of rainfall, G wi is twice as large as that in the absence of rainfall, which is illustrated in Figure 15a.

External Rainfall Disturbance Response
The transient responses of the GTDS to an unfavorable external rainfall disturbance are depicted in Figure 15. From the point of view that the occurrence of rainfall during the remediation will change the amount of liquid water diffusing into the heated soil firstly and affect the heating period and the energy consumption further, Gwi with step disturbances is used to simulate the rainfall events. Assuming that the rainfall occurs at each stage and lasts for 4 days, as well as that during periods of rainfall, Gwi is twice as large as that in the absence of rainfall, which is illustrated in Figure 15a. Under the action of the step disturbances in Gwi, the corresponding results are displayed in Figure 15b-d. As for the rainfall in Phase-1, owing to the fact that the amount of liquid water diffusion is extremely small, it has very little impact on vts and Gg, but leads to a small increase in vθs. When the rainfall occurs in Phase-2, vts remains 0 ℃/h with no change; vθs exhibits a positive 1.4% fluctuation for 12 h at the beginning of the rainfall and a negative 2.1% fluctuation for 14 h at the end of the rainfall; at the same time, Gg presents an increase from 0.00181 kg/s to 0.00185 kg/s in the early rainfall and a decrease from 0.00189 kg/s to 0.00183 kg/s with the end of the rainfall. Confronting 4 days of rainfall in Phase-3, vθs is fixed to 0 m 3 /(m 3 ·h) which is unaffected by the change of Gwi; vts exhibits a negative 11.6% fluctuation for 15 h when the rainfall begins and a positive 15.5% fluctuation for 18 h as the rainfall is over; meanwhile, Gg shows a rise from 0.00077 kg/s to 0.00086 kg/s with the occurrence of the rainfall and a decline from 0.00115 kg/s to 0.00106 kg/s with the end of the rainfall. Furthermore, the total mass of the natural gas is elevated by 82 kg under the action of rainfall.
In general, when confronting an unfavorable disturbance in the mass flow rate of liquid water diffusion caused by the external rainfall, the proposed control strategy demonstrates a high capacity Under the action of the step disturbances in G wi , the corresponding results are displayed in Figure 15b-d. As for the rainfall in Phase-1, owing to the fact that the amount of liquid water diffusion is extremely small, it has very little impact on v ts and G g , but leads to a small increase in v θs . When the rainfall occurs in Phase-2, v ts remains 0 • C/h with no change; v θs exhibits a positive 1.4% fluctuation for 12 h at the beginning of the rainfall and a negative 2.1% fluctuation for 14 h at the end of the rainfall; at the same time, G g presents an increase from 0.00181 kg/s to 0.00185 kg/s in the early rainfall and a decrease from 0.00189 kg/s to 0.00183 kg/s with the end of the rainfall. Confronting 4 days of rainfall in Phase-3, v θs is fixed to 0 m 3 /(m 3 ·h) which is unaffected by the change of G wi ; v ts exhibits a negative 11.6% fluctuation for 15 h when the rainfall begins and a positive 15.5% fluctuation for 18 h as the rainfall is over; meanwhile, G g shows a rise from 0.00077 kg/s to 0.00086 kg/s with the occurrence of the rainfall and a decline from 0.00115 kg/s to 0.00106 kg/s with the end of the rainfall. Furthermore, the total mass of the natural gas is elevated by 82 kg under the action of rainfall.
In general, when confronting an unfavorable disturbance in the mass flow rate of liquid water diffusion caused by the external rainfall, the proposed control strategy demonstrates a high capacity of Energies 2020, 13, 642 23 of 28 rapid response and self-adaptation. Additionally, the rainfall disturbance will result in the increase of the overall consumption of natural gas to some extent.

Conclusions
(1) For the purpose of actively adjusting the restoration period and reducing the energy consumption of polluted soil remediation, a natural gas thermal desorption system (GTDS) with a planned heating control strategy is proposed.
(2) The heating performance of the target site was significantly influenced by the natural gas flow. Increasing the mass flow rate of the natural gas lifted the stable temperature of the soil and shortened the settling time.
(3) The excess air coefficient exerted a large influence upon the soil's temperature and VWC characteristics. Smaller excess air coefficient which is greater than 1.05 may lead to a higher stable temperature and a shorter heating period of the heated block.
(4) With the planned heating control strategy, the heating process could be accomplished with constant temperature rise rate in Phase-1, fixed VWC change rate in Phase-2 and unvaried temperature rise rate in Phase-3 on schedule by adjusting the natural gas flow. Provided that the desired objective variables of each stage were set as 0.556 • C/h, −0.697 × 10 −3 m 3 /(m 3 ·h) and 1.18 • C/h, respectively, the overall heating time was shortened by 50.3% (36.5 days) and the total natural gas consumption was reduced by 24% (1487 kg) compared with the open-loop reference condition.
In conclusion, this planned heating controlling strategy applied for the GTDS can achieve the remediation period adjustment and energy consumption reduction, which are extremely meaningful for contaminated-soil restoration. The combination control of natural gas flow and excess air coefficient may be the development trend of natural gas thermal desorption systems in the future. It is expected that the results can serve as a theoretical reference for other thermal remediation systems for soil pollution.    Table A1. Geometric parameters of the GTDS.

Parameters Value (mm)
Burner [35] Inner diameter 200 Outer diameter 210 Length 1000 Thermal well [25] Inner diameter of inner pipe 760 Thickness of inner pipe 5 Length of inner pipe 4700 Inner diameter of outer pipe 1400 Thickness of outer pipe 5 Length of outer pipe 5000