The use of greenhouses in arid areas, such as in the Arabian Peninsula, faces serious obstacles due to the extremely harsh environment; daytime overheating caused by intensive solar irradiance and night-time unfavourable undercooling are common weather conditions in desert areas [1]. By solving greenhouse operational problems (cooling in summer days and heating in winter nights), farmers can grow protected crops throughout the year.

Evaporative cooling (wet pad–fan systems) and fuel burners are commonly used for cooling and heating greenhouse air, whereby cold humid air (in summer) and warm air (in winter) are exhausted from the greenhouse and lost to the surrounding (useless air) during ventilation processes. Moreover, heating and evaporative cooling systems are an additional burden of energy and water consumption in greenhouses [2]. To overcome these difficulties, a semi-closed greenhouse has been suggested, where the inside greenhouse air is circulated in a closed loop and its cooling or heating are obtained by using sustainable energy resources. Looking for low-cost, energy- and water-saving techniques for heating/cooling inside greenhouse air is an essential priority for sustainable agriculture

**Citation:** Al-Helal, I.; Alsadon, A.; Marey, S.; Ibrahim, A.; Shady, M.; Abdel-Ghany, A. Geothermal Energy Potential for Cooling/Heating Greenhouses in Hot Arid Regions. *Atmosphere* **2022**, *13*, 105. https:// doi.org/10.3390/atmos13010105

Academic Editors: Zengyun Hu, Xuguang Tang and Qinchuan Xin

Received: 10 December 2021 Accepted: 5 January 2022 Published: 10 January 2022

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

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

and the development of the rural areas in arid regions [3]. Therefore, a semi-closed greenhouse coupled with an earth–air heat exchanger (EAHE) was suggested in the Kingdom of Saudi Arabia (KSA). The EAHE is expected to be an effective way for heating/cooling greenhouses in the deserts characterised by hot arid climates. This is because the extracted low-grade geothermal energy (fluid temperature < 100 ◦C) as a sustainable energy resource is applicable for different purposes, such as desalination; drying technology; heating and cooling load production, using single-effect lithium bromide absorption chillers; agriculture greenhouses; and domestic use [4]. On the other hand, high-grade geothermal energy is applicable for power generation based on the Rankin cycle concept [5,6].

Even though the KSA is considered one of the most geothermally active countries in the Middle East [7], there is a lack of studies that address the issues of geothermal energy potential and its applications in the KSA [8]. Very few studies can be found that discuss the geothermal activity in the KSA from a geological science point of view [7,8]. Moreover, the geothermal energy potential, as well as the underground temperature profile in the shallow zone (0~8 m), has never been measured or evaluated in the KSA desert. Specifically, accurate information about the vertical distribution of soil temperature and geothermal energy potential in the shallow zone are still missing relatively to the KSA desert.

It is well known that the earth strongly absorbs solar energy during the daytime and stores a considerable amount of thermal energy at a particular depth. This is mainly attributable to the high solar irradiance intensity and large heat capacity of the soil. According to the diurnal variation in solar irradiance and ambient air temperature, the maximum temperature oscillation occurs at the Earth's surface, and it varies with the depth inside the Earth. Many researchers have measured the temperature distribution with depth and found that the oscillation of temperature below the ground attenuated with depth and became constant at a particular depth. For example, at a 0.18-m depth, the daily average soil temperature was 17.84 ◦C in the winter and 32.87 ◦C in the summer in New Delhi, India [9]. Moreover, the annual average temperatures of soil were measured at a 4-m depth and under different conditions of the ground surface in New Delhi, India. It was 29 ◦C when the ground surface was exposed to solar radiation, 19 ◦C when the ground surface was wetted, and 17 ◦C when the ground surface was wetted and shaded [9]. Different values of the annual average temperature of different soils in different places worldwide have been measured at different depths and reported in the literature, e.g., 23.45 ◦C at a 4-m depth in Las Vegas, USA [10]; 27–28 ◦C at a 2.5-m depth under warm humid weather conditions in Mexico City [11]; 18.7 ◦C at a 2-m depth and 20 ◦C at a 4-m depth in Brazil [12]; 25 ◦C at a 2-m depth in Bhopal, India [13]; 11.5 ◦C in winter and 17.5 ◦C in summer at a 2-m depth in Tianjin, China [14]; and 17.6 ◦C at a 3.6-m depth in Shouguang City, China [15]. In a typical arid climate (the desert in the south part of Algeria), the underground temperature was measured, in the hottest month in summer, to be 30 ◦C at a 2-m depth and 27 ◦C at a 5-m depth [16]. These temperatures are usually defined as ground undisturbed temperatures (GUTs). The GUT at a certain depth in the ground is nearly constant during the day and night and throughout the year. This depends on the groundwater level, the physical/chemical properties of soil and the ground surface conditions (mainly solar radiation flux and ambient air temperature). For greenhouse applications, the GUTs reported in the literature are reasonable for heating and cooling greenhouses; further, they can be used for indoor space heating and cooling (e.g., residential buildings, poultry houses, livestock houses, etc.) [17]. However, the cost of digging to the optimum depth for availing free geothermal energy should be considered; this depends on the type and nature of the soil. Besides the sustainability and low-cost and/or free geothermal energy for heating/cooling greenhouses, it can also be recognized that the thermal load levelling of the underground temperature at a depth from 2.5 m to 4 m is very low and is required for plants' healthy growth, irrespective of any climatic condition [17,18]. To use the geothermal energy potential for cooling/heating greenhouses in a closed loop, Polyvinyl Chloride (PVC) or High-density Polyethylene (HDPE) pipes are an optimum choice due to their low cost, low heat capacity and physical/chemical properties when they are buried at an optimum depth

in the ground. In greenhouse applications, buried pipes are to carry the exhausted humid air from the greenhouse; then, the air flowing through the pipes is either heated or cooled according to the requirement. This depends on the temperature difference between the inside greenhouse air and the inside surface temperature of the buried PVC or HDPE pipes. The mass flow rate, which depends on the cross-sectional area of the pipe and air velocity, can be optimized for a given greenhouse volume. This arrangement is generally referred to as the closed-loop earth–air heat exchanger (CL-EAHE). Such an arrangement has been implemented and evaluated for heating/cooling greenhouses by many researchers [17–23]; however, a survey of the literature revealed that there is a lack of and unclear information about CL-EAHEs operated in gravel–sand soil and arid climates (such as in the KSA). Moreover, because the buried pipes of CL-EAHEs carry humid greenhouse air, water vapour condensation inside the pipes is possible during the cooling process. This possibility needs to be examined for greenhouse applications. Another arrangement, called open-loop earth– air heat exchanger (OL-EAHE), is out of the scope of the present study. This is because, in the OL-EAHE, the preconditioned greenhouse air is usually discharged to the surrounding area outside the greenhouse (i.e., energy and water vapor losses) during the heating and cooling processes of inside greenhouse air. Accordingly, the main objectives of this study are (i) to experimentally measure the temperature profile at different depths inside the ground at the King Saud University campus (as a desert in an arid climate) to determine the optimum depth for burying EAHE pipes; (ii) to examine the possibility of condensation inside the buried pipes; and (iii) to evaluate the maximum heating and cooling potential of the geothermal energy in the KSA for possible use in greenhouses and other domestic and residential building applications.

#### **2. Methodology**

#### *2.1. Study Area*

The experiment was conducted at the Agricultural Research and Experiment Station, Agriculture Engineering Department, King Saud University (Riyadh, Saudi Arabia; 4643 E longitude and 2438 N latitude). The soil at the site of the experiment, as well as that of most of arid regions such as the Arabian Peninsula, is composed of gravelly sand and dries up to a 5-m depth or more. Therefore, mechanical digging was used to prepare a hole with a surface area of about 1 m2 and a depth of 3.5 m. In order to protect the cables of temperature sensors in the soil, these cables were collected to pass through a 5-cm diameter PVC pipe, installed vertically and fixed at the bottom of the hole with a concert block (Figure 1).

**Figure 1.** Schematic diagram showing the hole dug into the ground and the locations of thermocouples used to measure underground temperatures (*Tg*) and ambient temperature (*Tam*).

#### *2.2. Measuring Underground Temperature Profiles*

To determine the optimum depth at which the CL-EAHE pipes should be buried, the temperature profile of underground soil was measured at different depths. The optimum location for burying pipes is where the soil temperature became constant, unaffected by the diurnal and seasonal variations in the climatic parameters (GUT). The experiment was conducted during the period from 1 November 2020 to 31 October 2021, to measure the temperatures at the depths of 0.1, 0.5, 0.9, 1.3, 1.7, 2.1, 2.5 and 3 m, in addition to ambient air temperature. Digging up to only a 3.5-m depth was considered because more digging is quite hard and expensive due to the hard nature of soil in the desert of the KSA (i.e., mixture of gravelly sand and rocks). Moreover, the annual variation in soil temperature, *T*g (*T*g,max–*T*g,min) at the 3-m depth was expected to be low. Eight small holes were made vertically in the pipe wall (Figure 1) at the levels of measurement (3, 2.5, 2.1, 1.7, 1.3, 0.9, 0.5 and 0.1 m below the ground surface). These holes were to create a path for the allocation of the cables of the temperature sensors and to protect the cables to reach the data logger above the ground via the pipe in a safe (Figure 1). The sensors were inserted into the soil 0.8 m apart from the vertical pipe. The backfilling of the ground hole was performed in steps. The first filling was up to a 3-m depth; then, the temperature sensor was fixed carefully on the soil surface; after that, the second backfilling was made at the level of a 2.5-m depth; then, the temperature sensor was inserted, etc. Thermocouple sensors (wzp-035 Pt100/k; Shenzhen More-Suns Electronics Co., Ltd., Shenzhen, China) were used to measure the underground temperature at the specified locations. The thermocouple had a precision of ±0.1 ◦C and a temperature range of 0~85 ◦C. Measurements were taken every 5 min, averaged at every 15 min and recorded in a COMBILOG-1022 data logger (32 channels; Theoder Friedrichs & Co., Schenefeld, Germany).

Due to difficulties and to the high cost of digging, the vertical variation in *Tg* at different depths and at any time was calculated using the Kasuda formula [24]. Assuming homogeneous soil with a constant thermal diffusivity (αs), the monthly averaged soil temperature (*Tg*) at any depth (*z*) and month number (*t*) can be estimated by using the following formula:

$$T\_{\mathcal{S}}(z,t) = T\_{\text{m}} - T\_{\text{amp}} \times \exp\left[-z\sqrt{\left(\frac{\pi}{12a\_{\text{s}}}\right)}\right] \times \cos\left\{\frac{2\pi}{12}\left[t - t\_{\text{o}} - \frac{z}{2}\sqrt{\frac{12}{\pi a\_{\text{s}}}}\right]\right\} \tag{1}$$

where *Tg* (*z,t*) is the soil temperature at depth *z* and month number *t*; *T*m is the annual mean soil surface temperature (◦C); *T*amp is the amplitude of the soil surface temperature [(max − min)/2] in (◦C); *z* is the depth below the ground surface (m); αs is the ground soil thermal diffusivity (m2/month); *t* is time (the month number that the ground temperature is calculated for); *to* is the time shift (month of the year of the lowest ground surface temperature). For the desert of the KSA, the value of αs was estimated as 2.736 (m2/month) and thermal conductivity of soil as 2.2–2.8 (W m−1 ◦C−1) [7,8]. From 1 November 2020 to 31 October 2021, the values of *T*m and *T*amp were estimated (from measurements) to be 30 and 12 ◦C, respectively. Moreover, in January, the lowest ground surface temperatures were recorded; then, *to*, in Equation (1), is equal to one.

#### *2.3. Possibility of Condensation in the Buried CL-EAHE Pipes*

In the summer months, the inner surface temperature of the buried EAHE pipes is expected to be much lower than the temperature of the humid air exhausted from the greenhouse. This makes the geothermal effect for cooling inside greenhouse air possible. In summer, when humid air flowing through the buried pipes may cool to its dew-point temperature (*Tdp*) through contact with the inner surface of the pipe, which is cooler than the air, water vapour condenses on the pipe surface. In this case, design considerations should be taken into account to collect the condensed water from the buried EAHE pipes. This may increase the EAHE's cost; therefore, an accurate calculation for *Tdp* is important for the appropriate design of a low-cost CL-EAHE. To examine the possibility of condensation, an experiment was conducted to measure the dry-bulb temperature (*Td*) and relative humidity (*RH*) of the ventilated humid air at the greenhouse outlet (exactly before the exhaust fans). A crop-free greenhouse with a floor area of 165 m2 was used; it was covered with a double-layered polycarbonate sheet of an 8.15-mm thickness. The greenhouse was mechanically ventilated using two exhaust fans, each with an airflow rate of 350 m3 min−1, operated in summer (July 2021) and in winter (December 2020) to obtain data for two extreme weather conditions, in Riyadh, KSA. *Td* and *RH* were measured using a DMA033 Thermo-hygrometer (LSI-Lastem, Milano, Italy). The parameters were measured every 1 min, averaged at every 15 min and recorded in a data logger (CR23X Micro logger, Campbell Scientific, Inc., Oldenburg, Germany). The well-known approximate formulation used to calculate the dew-point temperature (*Tdp*) is based on Magnus' formula [25], for which the measured *Td* (in ◦C) and *RH* (%) are required; then, the saturated water vapour pressure (*Ps* in Pascal) corresponding to *Td* is given, by [25], as follows:

$$P\_s = 610.78 \times \exp\left(\frac{17.2694 \times T\_d}{(T\_d + 238.3)}\right) \tag{2}$$

The actual water vapour pressure (*Pa*, in Pascal) is given by the following:

$$P\_a = P\_s \times RH(\%) / 100\tag{3}$$

$$S = \ln(P\_a / 610.78) \tag{4}$$

$$T\_{dp} = \text{S} \times 238.3 / (17.294 - \text{S})\tag{5a}$$

Another simple approximation is used to calculate *Tdp* ( ◦C) for *RH* > 50% and ±1 ◦C error [26], in the following form:

$$T\_{dp} = T\_d - \left(\frac{100 - RH}{5}\right) \tag{5b}$$

These approximations are commonly used to calculate *Tdp*; however, in arid climates (where *Td* is very high and *RH* is very low), the validity of these approximations ((Equation (5a,b)) needs to be examined specifically for the arid climate.

For greater accuracy, the saturation water vapour pressure (*Ps*) has been modified and Equation (5a) enhanced, becoming known as an Arden Buck equation [27], by which the modified saturation water vapour pressure (*Psm*, in Pascal) corresponding to a dry bulb temperature of air (*Td*, in ◦C) is given by the following:

$$P\_{\rm sun} = a \times \exp\left\{ \left( b - \frac{T\_d}{D} \right) \left( \frac{T\_d}{c + D} \right) \right\} \tag{6}$$

$$T\_{dp} = \frac{c \times \ln\left(\frac{RH}{100} \times \frac{P\_{\text{sun}}}{a}\right)}{b - \ln\left(\frac{RH}{100} \times \frac{P\_{\text{sun}}}{a}\right)}\tag{7}$$

The constants in Equations (6) and (7) were provided with a maximum error of ≤0.05% and *Td* from 0 up to 50 ◦C as follows: *a* = 611.21 Pa, *b* = 17.368, *c* = 238.88 ◦C and *D* = 234.5 ◦C [27].

#### *2.4. Geothermal Cooling/Heating Capacity*

In an EAHE, the PVC or HDPE pipes are buried permanently under the ground at a specified location, where the ground temperature (*Tg*) is annually stable and nearly constant; therefore, the thermal conditions can be characterized as steady state and equilibrium. Moreover, the thickness of an EAHE pipe is very modest (e.g., 4.25 mm) compared to the pipe surface area or the pipe diameter (203.2 mm). Hence, for simplicity, we assumed that the thermal resistance of the pipe material was negligible, and the inner surface temperature

of the pipe (*Ts*) was almost uniform and constant along the axial direction of pipe and equal, in most cases, to the ground temperature (*Tg*) at the pipe location. The cooling and heating capacity of an EAHE ( . *Qc*/*h*) is defined as follows:

$$
\dot{Q}\_{c/h} = \dot{m}\_a \mathbb{C}\_p \left( T\_{\text{out}} - T\_{\text{in}} \right) \tag{8}
$$

where . *ma* is the mass flow rate of air exhausted from the greenhouse (ventilation rate at *Tex*) and flowing through the EAHE pipes (kg s−1); *Cp* is the specific heat capacity of air (J kg−1 ◦C−1); and *Tin* and *Tout* are the temperatures of the inlet and outlet air (◦C) for the EAHE pipes. *Tin* is assumed to be the temperature of air exhausted from the greenhouse (*Tin* = *Tex* in ◦C) and it can be obtained directly from measurements. For a very long pipe, *Tout* is assumed to be equal to the inner surface temperature of the EAHE pipes (i.e., *Tou*t = *Ts* = *Tg*). In this case, Equation (8) represents the maximum cooling/heating capacity that can be provided by an EAHE (*Qc/h,max*). In a similar manner, the cooling/heating potential (*Qc/h*), in Joules, over a specified period is defined, by [17,18], as follows:

$$Q\_{\mathbf{c}/\hbar} = \sum\_{time} \dot{m}\_a \mathbb{C}\_p (T\_{out} - T\_{in}) \Delta t \tag{9}$$

where Δ*t* is the time interval (s). In fact, the outlet air temperature of EAHE pipes (*Tout*) is lower than the ground temperature (*Tg*) or the inner surface temperature of the pipes (*Ts*); this depends on the effectiveness of the EAHE. Therefore, an expression for *Tout* in terms of underground temperature (*Tg*) and inlet temperature of air (*Tin*) is given, by [13], as follows:

$$T\_{out} = T\_{\mathcal{S}} \left\{ 1 - \exp\left( -\frac{A\_s h\_f}{\dot{m}\_a \mathbb{C}\_p} \right) \right\} + T\_{in} \exp\left( -\frac{A\_s h\_f}{\dot{m}\_a \mathbb{C}\_p} \right) \tag{10}$$

where *As* is the inner surface area of a single pipe in the EAHE; . *ma* and *Cp* are the mass flow rate and specific heat of air flowing in a single pipe (kg s−1 and J kg−1 ◦C−1); and *hf* is the convective heat transfer coefficient between the inner surface of the pipe and the flowing air (W m−2 ◦C−1). *hf* is a critical parameter in the EAHE design; in such a case, it is reasonable to assume that the air flow in the EAHE pipes is fully developed. To adapt an appropriate correlation to calculate *hf*, Bisoniya [13] examined eight Nusselt number (*Nu*) correlations and recommended the following correlation for turbulent flows in a pipe with a smooth internal surface [13]:

$$N\_{\rm u} = \frac{f/8(R\_c - 1000)P\_r}{1 + 12.7\sqrt{(f/8)}\left(P\_r^{2/3} - 1\right)}\tag{11}$$

The friction coefficient *f* for turbulent flow is given by the following:

$$f = \left(1.82 \log R\_{\varepsilon} - 1.64\right)^{-2}, \left(2300 \le R\_{\varepsilon} < 5 \times 10^{6}\right) \text{ and } \left(0.5 < P\_{r} < 10^{6}\right) \tag{12}$$

where the Reynolds number (*Re* = *dvρ*/*μ*) and Prandtl number (*Pr* = *μCp*/*k*) are calculated for the air flowing through a pipe; *d* is the inner diameter of the pipe; *v* is the air velocity in the pipe; *ρ* is the air density; *μ* is the air viscosity; and *k* is the thermal conductivity of air. The unit of each parameter was adjusted to give dimensionless numbers (*Nu*, *Re* and *Pr*).

By determining the value of *hf*, Equation (10) can be used as a design tool to optimize the number of pipes and the length and diameter of each pipe for the EAHE design.

#### **3. Results and Discussion**

*3.1. Underground Soil Temperature*

It is well known that the ground temperature is affected by the cyclic variation in climatic parameters such as solar radiation, ambient air temperature, wind speed, humidity, rainfall and snow, if any, etc. The interaction between the climatic parameters

and the ground is mainly at the ground surface; then, the cyclic variation in the ground temperature (*Tg*) decreases as the ground depth increases. At a certain depth, according to the chemical/physical properties of soil, *Tg* becomes nearly constant, unaffected by the diurnal and seasonal variations in the climatic parameters, and it is defined as the ground undisturbed temperature (GUT). The GUT is an essential parameter in designing an EAHE system. It is quite difficult to theoretically calculate the GUT correctly, because it depends on the soil parameters and the climatic conditions affecting the ground surface. Therefore, an experiment was conducted (Section 2) for measuring *Tg* at different ground depths during different seasons to determine the depth at which *Tg* became nearly constant throughout the year (GUT). To illustrate the daily cyclic variation in ambient air and soil temperatures during the period of the experiment (from 1 November 2020 to 31 October 2021), four days were selected to represent the extreme weather conditions, two days in cold winter (29–30 December) and two days in hot summer (30–31 July). The measured parameters (i.e., *Tam* and *Tg* at different depths) are depicted in Figure 2 (for winter days) and in Figure 3 (for summer days). In winter, the ground temperature (*Tg*) increased as the ground depth increased (Figure 2) and the opposite observation was noted in summer (Figure 3). In Figures 2 and 3, the cyclic variation in *Tg* is significant at a 0.1-m depth as affected by the diurnal variation in the climatic parameters (mainly the solar irradiance and *Tam*) in summer and winter. At a depth from ≥0.5 m up to 3 m, *Tg* became nearly constant in summer and in winter and was not affected by the diurnal variation in the climatic parameters. At a 3-m depth, *Tg* was about 29 ◦C in winter (Figure 2) and about 32 ◦C in summer (Figure 3) and a difference of 3 ◦C in the *Tg* value is acceptable between the hot summer and cold winter seasons. For more clarification and to find the optimum depth for burying the EAHE pipes, the vertical variation in the monthly average soil temperature (*Tg*) at different depths for the experimental site (Riyadh, KSA) is plotted in Figure 4 and the annual variation in temperatures of ambient air (*Tam*) and soil (*Tg*) at the 0.5-m and 3-m depths is plotted in Figure 5. As illustrated in Figures 4 and 5, at the 3-m depth, the maximum annual variation (max – min) in *Tg* was about 5 ◦C. The lowest value of *Tg* was around 27 ◦C in February and March and the highest value of *Tg* was around 32 ◦C in September and October. This annual variation in *Tg* (29–32 ◦C) is adequate for heating and cooling purposes and can be considered as the GUT and the depth of 3 m is adequate for burying the EAHE pipes. The annual variation in *Tg* is considerable at the 0.5-m depth; more stability of *Tg* was observed at the 3-m ground depth (Figure 5). According to Figure 5, the heating effect of a supposed EAHE in winter is expected to be relatively higher than its cooling effect in summer; this depends on the temperature difference between *Tam* and GUT (*Tg* at the 3-m depth).

**Figure 2.** Diurnal variation in ambient temperature (*Tam*) and underground soil temperature (*Tg*) at different depths (*z*) in winter season (29–30 December 2020).

**Figure 3.** Diurnal variation in ambient temperature (*Tam*) and underground soil temperature (*Tg*) at different depths (*z*) in summer season (30–31 July 2021).

**Figure 4.** Vertical variation in the monthly average ground temperature (*Tg*) at different depths (*z*) for Riyadh region, KSA, in the years 2020–2021.

As farther digging in the desert of arid regions is quite difficult and expensive, we used Equation (1) to predict the monthly average values of *Tg* at different depths. To calculate the values of *T*m and *T*amp in Equation (1), the measured *Tg* at the 0.1-m depth was considered as the soil surface temperature. The predicted results for 12 months and depths of >3 m are plotted in Figure 6. As illustrated in Figure 6, at a depth (*z*) greater than 13 m, *Tg* becomes constant at 30 ◦C throughout the year. As the arid climatic condition is almost similar every year in the deserts of KSA, the GUT can be considered constant at 30 ◦C throughout the year at a depth greater than 13 m.

**Figure 5.** Annual variation in temperatures of ambient air (*Tam*) and ground soil (*Tg*) at 0.5- and 3.0-m depths for Riyadh, Saudi Arabia.

**Figure 6.** Ground temperature variation (monthly average) with depth, predicted using Equation (1) for the KSA desert.

#### *3.2. Possibility of Condensation in the EAHE Pipes*

The proposed CL-EAHE was applied in a greenhouse for cooling/heating purposes in the summer and winter seasons. In the cooling/heating processes, humid air is exhausted from the greenhouse (at *Tex* and *RHex*); it flows in the buried pipes, then the condensation of water vapour in the pipes is possible during cooling. The occurrence of condensation requires extra design consideration to collect the condensed water from the pipes. The most accurate correlation (Equation (7)) was used to calculate the dew-point temperature (*Tdp*) of the exhausted air (having *Tex* and *RHex*) for the four previously selected days (hot and dry in summer and humid and cold in winter). No heating or cooling was applied to the greenhouse to represent the actual practical situations. Figures 7 and 8 illustrate the diurnal variation in the measured *Tex* and *RHex* for the exhausted air from the greenhouse, the measured GUT (*Tg* at the 3-m depth) and the calculated *Tdp* (using Equation (7)). A considerable difference could be observed between *Tex* and *Tdp* in summer (Figure 8) and in winter (Figure 7). For winter heating (Figure 7), *Tex* increased, as affected by the warm inner surface of the EAHE pipes; if the upper limit of *Tex* was to reach the GUT, then the gap

between the *Tex* and *Tdp* of the air would increase and the condensation of water vapour would be impossible. For summer cooling (Figure 8), if *Tex* was to decrease to its lower limit GUT (i.e., *Tg* at the 3-m depth), there would still be a large difference between *Tex* and *Tdp* (Figure 8); therefore, condensation would never take place or would be impossible. Moreover, a large number of combinations for *Tex* and *RHex* was used to calculate *Tdp* using Equation (7) and the resulting values of *Tdp* do not exceed 20 ◦C; in each combination, a large difference between *Tex* and *Tdp* remained. Accordingly, condensation in CL-EAHE pipes would be impossible throughout the year for the pipes buried at a 3-m depth in the Riyadh area, KSA.

**Figure 7.** Diurnal variation in estimated dew-point temperature (*Tdp*) of air ventilated from a greenhouse at *Tex* and *RHex* and flowed through EAHE pipes buried at the GUT in winter season (29–30 December 2020).

**Figure 8.** Diurnal variation in estimated dew-point temperature (*Tdp*) of air ventilated from a greenhouse at *Tex* and *RHex* and flowed through EAHE pipes buried at the GUT in summer season (30–31 July 2021).

Similarly, *Tdp* was calculated using Equation (5a,b) to check the validity of these approximations to be used to calculate *Tdp* for an arid climate. The percentage error, *ER* (%), was calculated in each case as *ER* (%) = abs{(Taccurate-Tapproximate)/Taccurate} × 100. The values of *ER* (%) were calculated for a wide range of *Td* and *RH* and are plotted in Figure 9a,b for Equation (5a) and in Figure 10 for Equation (5b). Based on Figure 9a,b, Equation (5a) can be used to calculate *Tdp* with an error less than 1% for *RH* > 40% and *Td* < 30 ◦C. However, Equation (5b) can be used only for *RH* > 60% with an error of about 3–4%; for low *RH* values (as in the arid climate), such an approximation (Equation (5b)) cannot be used. Accordingly, for arid climatic conditions (as in the Arabian Peninsula region), such approximations (Equation (5a,b)) are not recommended and Equation (7) is the appropriate correlation to determine *Tdp* correctly.

**Figure 9.** Percentage error, *ER* (%), in the approximated value of *Tdp* (estimated by using Equation (5a)) as affected by (**a**) relative humidity (*RH*) and (**b**) dry bulb temperature (*Td*).

**Figure 10.** Percentage error, *ER* (%), in the approximated value of *Tdp* (estimated by using Equation (5b)) as affected by relative humidity (*RH*).

#### *3.3. Maximum Cooling/Heating Capacity*

An EAHE's geothermal energy potential, as well as its cooling/heating capacity, depends mainly on the climatic parameters and the underground soil conditions. For EAHE pipes buried at a specified optimum depth, *Tg* equals the GUT. Under the steadystate thermal condition, the inner surface temperature of the pipe (*Ts)* is assumed to be uniform in the axial flow direction and by assuming the thermal resistance of the pipe material is negligible, the equality *Ts* = *Tg =* GUT can be assumed with insignificant error. If the pipe is long enough, the outlet air temperature (*Tout*) from the pipe can be assumed to be equal to *Ts* (*Tout* = *Ts* = *Tg* = GUT). In this case, the EAHE provides its maximum possible cooling/heating capacity. Therefore, in Equation (8), *Tin* and *Tout* were taken as *Tex* (exit air from the greenhouse) and GUT, respectively; . *Qc*/*h*,*max* was estimated per cubic meter of airflow. For winter cooling, a considerable amount of heat is expected to be added to the flowing air at around midnight (Figure 11). In addition, at around noon, in winter, operating the EAHE is not necessary and the transmitted solar radiation into the greenhouse is enough for warming up the inside air. In summer, a considerable amount of heat is expected to be removed from the flowing air at around noon; however, at around midnight, it is not necessary to operate the EAHE (Figure 11). Under the presumed ideal conditions and to estimate the maximum possible cooling/heating potential that the ground (at a 3-m depth) can provide in the cold winter and hot summer in the Riyadh region, the results are integrated in Figure 11. The expected maximum possible cooling and heating potential are 890 and 1000 MJ m−3 day−1, respectively. These values are promising to use EAHEs for different applications in the KSA and the Arabian Peninsula regions.

**Figure 11.** Diurnal variation in maximum cooling/heating capacity ( . *Qc*/*h*,*max*) that an EAHE can provide and the maximum cooling/heating potential per day.

#### **4. Conclusions**

This study is an attempt to provide critical information for EAHE designers and to evaluate the geothermal energy potential for possible cooling/heating greenhouses and indoor spaces in hot arid regions, for which two experiments were conducted. The main conclusion is summarized below.

The daily cyclic variation in the climatic parameters affected the ground temperature (*Tg*) up to a 0.5-m depth, after which *Tg* was constant daily; it increased with depth in winter and decreased with depth in summer.

At a 3-m depth, the annual variation in *Tg* was minor and it can be considered as the ground undisturbed temperature, GUT. The GUT value was 29 ◦C in winter and 32 ◦C in summer; these values are adequate for cooling/heating purposes in a climate in which the ambient air temperature drops below 10 ◦C on winter nights and exceeds 47 ◦C on summer days. For an EAHE operating under ideal conditions, the geothermal energy level can provide a maximum cooling/heating potential of 890/1000 MJ per m3 of flowing air per day.

In summer seasons, during the cooling process of greenhouse air in EAHE pipes distributed at a 3-m depth, the condensation of water vapor in the pipes could never take place.

In arid climates, approximate correlations (Equation (5a,b)) are not recommended to calculate the dew-point temperature of air (*Tdp*). However, Equation (7) is the appropriate correlation to determine *Tdp* correctly.

The geothermal energy potential is promising for cooling/heating applications in arid regions such as the Arabian Peninsula for sustainable development and environmental protection. Further research should be conducted to design, construct and operate an EAHE connected to a greenhouse and evaluate the system performance under different operating and climate conditions.

**Author Contributions:** Conceptualization and methodology, A.A.-G. and I.A.-H.; validation, A.A., S.M. and A.I.; formal analysis, A.A.-G., S.M. and I.A.-H.; investigation, A.A.-G.; resources, M.S.; data creation, M.S., S.M. and A.I.; writing—original draft preparation, I.A.-H., S.M. and A.I.; writing review and editing, A.A.-G.; visualization, M.S., S.M. and A.I.; supervision, A.A.-G.; project administration, I.A.-H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research study received no external funding.

**Acknowledgments:** The authors extend their appreciation to the Deputyship for Research and Innovation, Ministry of Education in Saudi Arabia, for funding this research work through the project number DRI-KSU-542.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**

#### **References**

### *Article* **Decomposition and Decoupling Analysis between HDI and Carbon Emissions**

**Dongju Li 1, Tongtong Shen 2,\*, Xi Wei 1 and Jie Li 1**

**Abstract:** The concept of low carbon is extended to the welfare dimension by considering the relationship between carbon emissions and the Human Development Index (HDI). This paper examines the decoupling between carbon emissions per capita and HDI and the welfare output of carbon emissions by using the data from 189 countries, from 1990 to 2019, as well as decomposes the drivers of the decoupling index and carbon emissions performance (CEP) in the example countries. The results show that most countries that achieve strong decoupling have very high human development, while the worst case is that a few countries with an extremely low human development achieved strong decoupling. Moreover, the status of strong decoupling in most countries is not stable, and there is a risk of transformation to another decoupling status. Although the CEP of most countries has gradually improved, very few countries have high CEP. Economic development effects are the primary inhibitor to achieving and maintaining strong decoupling in example countries. The main drivers of CEP improvement are the carbon productivity effects in the Czech Republic, Germany, and the United Kingdom, and the economic development effects in South Korea and Turkey. The main factors inhibiting the increase of CEP are the energy intensity effect in the Czech Republic, Germany, and the UK, and the welfare effect in South Korea and Turkey. These effects are all related to GDP. Economic activity broadly affects the decoupling index and CEP. Recommendations for maintaining HDI growth and reducing carbon emissions are made for countries with different human development.

**Keywords:** HDI; decoupling index; carbon emission performance; LMDI

#### **1. Introduction**

In the context of global warming, the low-carbon economy with low energy consumption, low pollution, and low emissions have become the focus of global attention. In compliance with the *Kyoto Protocol* and the *Paris Agreement*, governments have taken actions to reduce carbon emissions, improve energy efficiency, and enhance their ability to cope with climate change. Despite this, global energy demand and CO2 emissions continue to increase, except in 2020 [1]. Ideally, we would like to see economic growth that is not dependent on the increase in carbon emissions, with decoupling between economic growth and carbon emissions. Therefore, revealing the decoupling process is conducive to finding critical breakthroughs in the path of global low-carbon economic development [2].

Low-carbon development is being achieved at the economic level, but it should also be extended to the welfare dimension. Based on Sen's capability approach, the Human Development Index (HDI) is constructed to measure the combined performance of a country in health, education, and income [3]. The increase in HDI is taken as a starting point for human welfare. To achieve higher human development is the pursuit of each country. However, realizing the right to human development requires the guarantee of energy consumption and the demand for carbon emissions [4]. The 30th anniversary

**Citation:** Li, D.; Shen, T.; Wei, X.; Li, J. Decomposition and Decoupling Analysis between HDI and Carbon Emissions. *Atmosphere* **2022**, *13*, 584. https://doi.org/10.3390/ atmos13040584

Academic Editors: Zengyun Hu, Xuguang Tang and Qinchuan Xin

Received: 25 February 2022 Accepted: 2 April 2022 Published: 5 April 2022

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

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

Human Development Report introduced a new index—Planetary pressures-adjusted HDI (PHDI), which emphasizes how human activities significantly influence the environment and ecology and place great stress on the planet. These pressures will react to human development in the form of extreme climate or species extinction, resulting in the decline or loss of HDI [5]. Therefore, studying the relationship between CO2 emissions and human development will provide policymakers with more evidence for achieving sustainable development [6].

This paper aims to analyze the decoupling between carbon emissions and HDI and the carbon welfare output of 189 countries from 1990 to 2019 and understand the drivers. First, the overall decoupling performance of 189 countries from 1990 to 2019 is analyzed. Furthermore, the evolution of decoupling trends is shown over six periods. However, the decoupling index and decoupling status are not comprehensive enough to explain the coupling relationship between carbon emissions and HDI. Therefore, the total carbon emissions and carbon emissions per capita were combined to analyze the welfare output of carbon emissions in 189 countries, namely carbon emission performance (CEP). Finally, several countries that maintained continuous decoupling and expansive negative decoupling over six periods are identified in the decoupling analysis. And LMDI analyzed the driving factors of the decoupling index and CEP.

This paper consists of five sections and is arranged as follows. The second section reviews and summarizes the relevant literature. The appropriate methods and data are described in Section 3. Section 4 shows the results of a series of methods. Section 5 summarizes the main conclusions and makes recommendations on maintaining a stable decoupling status.

#### **2. Literature Review**

#### *2.1. Review of Decoupling Model*

The decoupling model is widely used to analyze the coupling relationship between economy and resource consumption or environmental pollution. The two main decoupling models are the OECD decoupling model and the Tapio decoupling model. OECD analyzed the decoupling status of 30 member countries and confirmed the widespread decoupling between environment and economy in OECD countries [7]. Tapio adopted the concept of elasticity to dynamically reflect the decoupling relationship between variables and determined eight reasonable decoupling statuses, overcoming the difficulties in the base period selection of the OECD decoupling model [8]. Lin et al. used the OECD and Tapio decoupling model to evaluate the decoupling status between CO2 emissions and GDP in South Africa during 1990–2012 and found that the decoupling effect was gradually increasing [9]. The results of the two decoupling models have a good correlation, but the OECD decoupling model has limitations for analyzing the decoupling of sub-periods. To compare the decoupling trend between economic growth and carbon dioxide emissions in typical developed and developing countries from 1965 to 2015, Wu et al. used three decoupling models—the OECD decoupling model, the Tapio decoupling model, and the IGTX decoupling model [10]. The research has shown that the decoupling effect is stronger in developed countries than in developing countries, and the accuracy of the Tapio decoupling model is not easily limited by the length of the research period among the three decoupling models.

In the decoupling analysis between carbon dioxide and economy, the Tapio model and LMDI decomposition method are generally combined to discuss the driving factors affecting the decoupling status. Chen et al. embedded the LMDI decomposition formula in the Tapio decoupling model to study the driving factors affecting CO2 emissions in OECD countries from 2001 to 2015 [11]. The results showed that energy intensity and GDP per capita are the main factors affecting CO2 emissions, and the impact of technological factors is greater than that of non-technological factors. Yang et al. combined the LMDI decomposition model and Tapio decoupling model to analyze the determinants of decoupling between global carbon emissions and global economic growth during 2000–2017 [2]. It

is found that energy conservation, production efficiency, and energy structure optimization played a positive role in the global decoupling process, while GDP per capita and population expansion played a restraining role. In general, the application of the Tapio model and LMDI decomposition model is wide. Therefore, the Tapio model and LMDI decomposition model are used in this paper to identify each country's decoupling status and driving factors.

#### *2.2. Review of Carbon Emissions and HDI*

Most literature focuses on the relationship between energy consumption/environmental pollution and economic growth but ignores the relationship between energy/environment and human welfare. Economic growth is only one of the goals of human development. In essence, promoting human development is the ultimate goal. In the context of limited global energy resources, ensuring the effective use of energy and promoting human development is particularly important. Martínez and Ebenhack found a significant correlation between HDI and energy consumption per capita in 120 countries [12]. For some energy-poor countries with a low HDI, a small amount of energy assistance can make great strides in human development. Akizu-Gardoki et al. conducted a decoupling analysis of 126 countries from 2000 to 2014, based on Total Primary Energy Footprint and HDI, and the results showed that human welfare could be increased by reducing energy consumption [13]. According to the analytical framework of HDI, there is also research on how to formulate and implement global emission reduction policies. Pan believes that countries with a low HDI should be given priority in international negotiations on global climate change mitigation to guarantee their rights in achieving human development [4]. Hu proposed to implement global emission reduction policies according to the HDI of countries/regions and the principle that major carbon emitters take the lead in emission reduction [14]. Costa et al. used the positive correlation between HDI and CO2 emissions per capita to evaluate the cumulative CO2 emissions required to achieve high human development in developing countries and determine their respective emission reduction according to the HDI [15].

The decoupling relationship between carbon emissions and HDI is mainly studied at the levels of nation and district. Based on the Tapio decoupling model and LMDI decomposition, Chen et al. discussed the decoupling status between carbon emissions per capita and the HDI of four provinces in southwest China from 2000 to 2015 [16]. The results showed that the overall decoupling effect of southwest China was gradually increasing. Hossain and Chen demonstrated weak and strong decoupling between HDI and CO2 emissions in Bangladesh from 1990 to 2018, with economic activity being the primary driver of CO2 emissions, compared to other influencing factors [6]. Therefore, this paper expands the research perspective to the global dimension instead of studying a specific country or region. The decoupling relationship between carbon emissions and HDI can be studied from the perspective of multiple countries, and comparative analysis can be conducted among countries with the same data source.

#### *2.3. Review of Carbon Emission Performance*

Carbon productivity originated from the resource productivity in ecological economics and was later measured by consultants McKinsey & Company, in terms of GDP output per unit of CO2 emissions [17]. According to the definition of carbon productivity, carbon emission performance (CEP) measures the economic and social welfare output per unit of carbon emission [18,19]. Therefore, the welfare output of carbon emissions should be considered from the perspective of human well-being to broaden the research horizon of low-carbon development. Zhu and Liu measured the CEP of G20 countries, expressed by the ratio of HDI to carbon emission index—the average of the proportion of total carbon emission and carbon emission per capita after standardization [18]. Liu and Zhu and Hu et al., respectively, measured the carbon emission performance of G20 countries and 30 Provinces in China, but the carbon emission index only included carbon emissions per capita, without the important information of total carbon emissions [19,20]. The PHDI in the Human Development Report 2020 adjusts HDI through carbon emissions per capita and material footprint per capita to explain the pressure on the earth caused by human activities [5]. In this paper, the calculation method of Zhu and Liu and PHDI was used to construct the carbon emission performance index.

#### **3. Methods and Data**

#### *3.1. Analysis of Decoupling*

In economics, the decoupling model is usually used to study the coupling relationship between economic growth and resource consumption or environmental pollution, reflecting the sensitivity of resources and environmental pressure changes to economic changes [21]. This paper uses the Tapio decoupling model to analyze the decoupling relationship between carbon emissions per capita and HDI. Since HDI is a per capita indicator, carbon emissions per capita are selected instead of total carbon emissions. The equation is expressed as follows.

$$D(PC, HDI) = \frac{\% \Delta PC}{\% \Delta HDI} = \frac{\Delta PC/PC\_0}{\Delta HDI/HDI\_0} \tag{1}$$

In Equation (1), *D*(*PC*, *HDI*) represents the decoupling index between carbon emissions per capita (*PC*) and the Human Development Index (*HDI*), and %Δ*PC* and %Δ*HDI* represent the growth rate of carbon emissions per capita and HDI, respectively. Both Δ*PC* and Δ*HDI* represent the change of carbon emissions per capita and HDI from the base period to the end period respectively. Lastly, *PC*0 and *HDI*0 represent carbon emissions per capita and HDI at the base period, respectively.

According to the decoupling index, the decoupling status is divided into eight statuses: strong decoupling, weak decoupling, expansive coupling, expansive negative decoupling, weak negative decoupling, strong negative decoupling, recessive coupling, and recessive decoupling. The meanings of the eight decoupling statuses are detailed in the other literature [22].

#### *3.2. Construction of Carbon Emission Performance (CEP)*

According to the definition of CEP, HDI is used to represent welfare. The HDI data provided by the Human Development Data Center is non-continuous, and direct use may cause errors. According to the calculation method, the HDI of 189 countries from 1990 to 2019 was recalculated, but the calculation steps were omitted. Carbon emissions are expressed by Carbon Emission Index (CEI). CEI needs to consider two factors: total carbon emissions and carbon emissions per capita. Zhu and Liu standardized the ratio of total carbon emissions to the total global emissions as the total Carbon Emission Index [18]. However, the total global emissions are determined by the total carbon emissions of the sample countries, and the proportion of total carbon emissions is not stable. In this paper, the total carbon emissions and carbon emissions per capita are directly standardized. Construction of carbon emission performance is shown in Equation (2):

$$CEP = \frac{HDI}{CEI} \tag{2}$$

$$CEI = \frac{1}{2} \left[ \frac{Ln(\text{C}) - Ln(0)}{Ln(10, 490, 000) - Ln(0)} + \frac{Ln(PC) - Ln(0)}{Ln(68, 720) - Ln(0)} \right] \tag{3}$$

In Equation (3), *C* is carbon emission and *PC* is carbon emission per capita. The standardized method for calculating PHDI is used [5]. Zero is set as the minimum value. The maximum value is the highest value in all countries from 1990 to 2019. Therefore, the minimum value of total carbon emission and carbon emission per capita is 0. The maximum value of total carbon emission is 10.49 billion tons, observed in China in 2019, and the maximum value of carbon emission per capita is 68.72 tons, observed in Qatar

in 1997. It should be noted that logarithmic standardization is meaningful only if the value of total carbon emission and carbon emission per capita is greater than 1. Therefore, the total carbon emission is expressed in 10,000 tons, and the carbon emission per capita is expressed in kilograms.

#### *3.3. Analysis of Decomposition*

Logarithmic Mean Divisia Index (LMDI) decomposes the variables into several effects and analyzes the contribution of each effect to the change in the variable [23]. This method can show the decomposition results as additivity, complete decomposition, and no residual error. The LMDI method has been widely used to decompose drivers of carbon emission due to its variety of forms and simplicity of calculation [24,25]. This paper decomposes the decoupling index and CEP by combining the Kaya Identity and the LMDI method.

#### 3.3.1. Decomposition of Decoupling Index

Based on the idea of existing studies, the decomposition process of drivers affecting CO2 emission is determined as follows [6,16]:

$$PC = \frac{C}{F} \times \frac{F}{E} \times \frac{E}{GDP} \times \frac{GDP}{P} \tag{4}$$

where, *CF* = *C F* , *FE* = *F E* , *EG* = *E GDP* , *GP* = *GDP P*

In Equation (4), *C* represents total carbon emissions, *F* represents fossil fuel consumption (oil, natural gas, and coal), *E* represents primary energy consumption, *GDP* represents total economic output, and *P* represents the size of the population. *CF* represents the carbon emission intensity of fossil fuel. *FE* represents the structure of energy consumption. *EG* represents energy intensity. *GP* represents GDP per capita. Equation (4) can be transformed into:

$$PC = CF \times FE \times EG \times GP \tag{5}$$

According to the LMDI model, the change of carbon emissions per capita from the base period (0 period) to the reporting period (T period) is shown in Equation (6).

$$\begin{split} \Delta P \mathbf{C} = P \mathbf{C}\_t - P \mathbf{C}\_0 &= \mathbf{C} F\_t \times F E\_t \times E \mathbf{G}\_t \times G \mathbf{P}\_t - \mathbf{C} F\_0 \times F E\_0 \times E \mathbf{G}\_0 \times G \mathbf{P}\_0 \\ &= \Delta \mathbf{C} F + \Delta F E + \Delta E \mathbf{G} + \Delta G P \end{split} \tag{6}$$

$$
\Delta \text{CF} = \text{Ln} \left( \frac{\text{CF}\_t}{\text{CF}\_o} \right) \times \frac{\text{PC}\_t - \text{PC}\_0}{\text{Ln} (\text{PC}\_t) - \text{In} (\text{PC}\_0)} \tag{7}
$$

$$
\Delta FE = Ln \left(\frac{FE\_t}{FE\_0}\right) \times \frac{PC\_t - PC\_0}{Ln(PC\_t) - ln(PC\_0)} \tag{8}
$$

$$
\Delta EG = Ln \left(\frac{EG\_t}{EG\_0}\right) \times \frac{PC\_t - PC\_0}{Ln(PC\_t) - ln(PC\_0)} \tag{9}
$$

$$
\Delta GP = Ln \left(\frac{GP\_t}{GP\_0}\right) \times \frac{PC\_t - PC\_0}{Ln(PC\_t) - ln(PC\_0)} \tag{10}
$$

In Equation (6), Δ*CF* represents the carbon intensity of fossil fuel effect, Δ*FE* represents the energy structure effect, Δ*EG* represents energy intensity of fossil fuel effect, and Δ*GP* represents economic development effect. The four effects from the base period to T period are calculated based in Equations (7)–(10).

By combining the Tapio decoupling model with the LMDI method, the decoupling index from the base period (0 period) to the reporting period (T period) can be expressed by Equation (11).

$$D(PC, HDI) = \frac{\% \Delta P\%}{\% \Delta HDI} = \frac{\Delta P\% / P\%\_0}{\Delta HDI / HDI\_0} = \Delta P\% \times \frac{HDI\_0}{\Delta HDI \times P\%\_0}$$

$$= (\Delta CF + \Delta FE + \Delta EG + \Delta GP) \times \frac{HDI\_0}{\Delta HDI \times P\%\_0}$$

$$\begin{aligned} \Delta = \Delta \text{CF} \times \frac{HDI\_0}{\Delta HDI \times PC\_0} + \Delta FE \times \frac{HDI\_0}{\Delta HDI \times PC\_0} + \Delta EG \times \frac{HDI\_0}{\Delta HDI \times PC\_0} + \Delta GP \times \frac{HDI\_0}{\Delta HDI \times PC\_0} \\ &= D\_{CF} + D\_{FE} + D\_{EG} + D\_{GP} \end{aligned} \tag{11}$$

In Equation (11), the decoupling index of carbon emissions per capita and HDI is divided into four components. *DCF*, *DFE*, *DEG*, and *DGP* represent the decoupling index of Δ*CF*, Δ*FE*, Δ*EG*, and Δ*GP*, respectively.

#### 3.3.2. Decomposition of CEP

Similarly, the Kaya Identity and the LMDI method are used to decompose the driving factors of CEP. The decomposition form is shown in Equation (12).

$$CEP = \frac{CEP}{GDP} \times \frac{GDP}{C} \times \frac{C}{F} \times \frac{F}{E} \times \frac{E}{GDP} \times \frac{GDP}{P} \times P \tag{12}$$

where, *CG* = *CEP GDP* , *GC* = *GDP C* , *CF* = *C F* , *FE* = *F E* , *EG* = *E GDP* , *GP* = *GDP P* , *P* = *POP*.

In Equation (12), *CG* represents the contribution to CEP per unit of economic growth. *GC* represents the economic output per unit of carbon emissions, i.e., carbon productivity. *CF*, *FE*, *EG*, and *GP* are the same as in Equation (4), respectively. *P* denotes the size of the population. Equation (12) is simplified into the following form:

$$
\mathcal{C}EP = \mathcal{C}G \times \mathcal{G}\mathcal{C} \times \mathcal{C}F \times FE \times EG \times GP \times P \tag{13}
$$

The change of CEP from 0 period to T period is expressed as Equation (14).

$$\Delta \mathbb{C}EP = \mathbb{C}EP\_l - \mathbb{C}EP\_0 = \Delta \mathbb{C}G + \Delta \mathbb{C}G + \Delta \mathbb{C}F + \Delta FE + \Delta EG + \Delta GP + \Delta P \tag{14}$$

$$\Delta \text{CG} = \text{Ln}\left(\frac{\text{CG}\_t}{\text{CG}\_o}\right) \times \frac{\text{CEP}\_t - \text{CEP}\_0}{\text{Ln}\left(\text{CEP}\_t\right) - \ln\left(\text{CEP}\_0\right)}\tag{15}$$

$$
\Delta GC = Ln \left(\frac{GC\_t}{GC\_0}\right) \times \frac{CEP\_t - CEP\_0}{Ln(CEP\_t) - ln(CEP\_0)} \tag{16}
$$

$$
\Delta CF = Ln \left(\frac{CF\_t}{CF\_0}\right) \times \frac{CEP\_t - CEP\_0}{Ln(CEP\_t) - ln(CEP\_0)} \tag{17}
$$

$$
\Delta FE = Ln \left(\frac{FE\_t}{FE\_o}\right) \times \frac{CEP\_t - CEP\_0}{Ln(CEP\_t) - ln(CEP\_0)} \tag{18}
$$

$$
\Delta EG = Ln \left(\frac{EG\_t}{EG\_0}\right) \times \frac{CEP\_t - CEP\_0}{Ln(CEP\_t) - ln(CEP\_0)} \tag{19}
$$

$$
\Delta GP = Ln\left(\frac{GP\_t}{GP\_0}\right) \times \frac{CEP\_t - CEP\_0}{Ln(CEP\_t) - ln(CEP\_0)}\tag{20}
$$

$$\Delta P = Ln \left(\frac{P\_t}{P\_0}\right) \times \frac{CEP\_t - CEP\_0}{Ln(CEP\_t) - ln(CEP\_0)} \tag{21}$$

In Equation (14), Δ*CG* represents the welfare effect, Δ*GC* represents the carbon productivity effect, and Δ*P* represents the population effect. The terms Δ*CF*, Δ*FE*, Δ*EG*, and Δ*GP* are the same as in Equation (6), respectively. The seven effects from the base period to T period are calculated based on Equations (15)–(21).

#### *3.4. Data*

Frequently used data on carbon emissions are taken from BP statistics.  The carbon emissions data used in this paper are from the Global Carbon Project, which mainly refers to the carbon emissions generated by burning fossil fuels within the national territory [26]. UNDP also uses carbon emissions data from the Global Carbon Project

to calculate PHDI [5]. *BP Statistical Review of World Energy 2021* provides a long period of energy consumption data for the decomposition analysis of countries with certain characteristics [1]. The GDP data derived from World Development Indicators is updated in March 2021 [27]. Population data were derived from the World Population Prospects: The 2019 Revision [28]. Raw data on health, education, and income were obtained from the Human Development Data Centre to recalculate the HDI [29]. The variable definitions and sources are shown in Table 1.

**Table 1.** Variable definition and source.

#### **4. Result**

#### *4.1. Relationship between Carbon Emissions Per Capita and HDI*

Figure 1 depicts the corresponding relationship between CO2 emissions per capita and the HDI for four years. There is a strong correlation between CO2 emissions per capita and the HDI. The figures show a power–exponential relationship, that is, with the increase of HDI, the increase of carbon emissions per capita is gradually accelerating. Carbon emissions per capita of countries with extremely low HDI are minor, and the growth is flat. When countries achieve medium, high, and very high human development, the growth rate of carbon emissions per capita begins to accelerate. This indicates that the demand for carbon emissions and the carbon emissions per capita will gradually increase for the country to achieve a higher level of human development. The non-linear relationship between CO2 emissions per capita and HDI is consistent with Costa et al., who found an exponential relationship between CO2 emissions from fossil fuels and HDI [15].

It can also be seen from Figure 1 that the CO2 emissions per capita of countries with extremely low HDI are relatively close, with little difference between countries. As the HDI of countries gradually increases, differences in carbon emissions per capita become apparent. The differences in emissions per capita among countries are most remarkable at very high human development. Since 1990, emissions per capita by the poorest half of the world's population have only slightly increased, from 1.2 tons to 1.6 tons [30]. For countries with high levels of human development, their carbon emissions per capita are different, and they follow different development paths. Some countries will achieve high human development at the cost of large carbon emissions. Some countries will aim at low-carbon development, improving emission-reduction technologies, and developing new energy sources to achieve high human development.

**Figure 1.** Relationship between HDI and CO2 per capita emissions. Source: Authors' own construction. Note: The *p*-values for the four exponential functions are all less than 0.01 and pass the significance test. (**a**) 1990; (**b**) 2000; (**c**) 2010; (**d**) 2019.

#### *4.2. Analysis of Decoupling between Carbon Emissions Per Capita and HDI*

#### 4.2.1. Decoupling between Carbon Emissions Per Capita and HDI, 1990–2019

It is the right of every country to achieve high human development. Carbon emissions should be targeted to increase the HDI. If carbon emissions per capita increase without increasing the HDI, this may indicate that no effort has been made to reduce carbon emissions. Countries with strong decoupling and very high human development are considered to have achieved positive decoupling. However, countries that do not achieve very high HDI and strong decoupling are referred to as negative decoupling.

As shown in Figure 2, no country experienced negative growth HDI from 1990 to 2019. HDI data for 1990 are missing for 46 countries, and 59 countries out of the 143 are strongly decoupling. In addition, 40 countries are positively decoupling whereas 19 countries are negative decoupling. Negative decoupling is not a good thing for countries with extremely low HDI. These countries are the Central African Republic, the Democratic Republic of Congo, and Yemen. Their HDI is still at low human development after 29 years, due to constant war and political unrest. Because of their weak industrial base and poor infrastructure, they are strongly decoupling with low carbon emissions and low energy consumption. Much investment is needed to develop in the future due to the low starting point.

**Figure 2.** Decoupling between carbon emissions per capita and HDI, 1990–2019. Source: Authors' own construction using Excel.

In total, 63 countries are expansive negative decoupling, meaning carbon emissions per capita are growing much faster than HDI. Most of the 63 countries are at high or medium human development, and 11 countries with very high HDI are expansive negative decoupling, such as South Korea, Kuwait, and Qatar. Achieving high HDI growth comes at the cost of high carbon emissions per capita.

#### 4.2.2. Evolution of the Decoupling between Carbon Emissions Per Capita and HDI

The whole observation period from 1990 to 2019 is divided into six sub-periods. The decoupling index between carbon emissions per capita and HDI is calculated to analyze the change in the decoupling status in each period. As shown in Figure 3, the decoupling status of countries changed significantly over the six periods. Most countries have maintained positive HDI growth. In general, the main decoupling statuses are strong decoupling, weak decoupling, expansive coupling, and expansive negative decoupling. The proportion of countries with these four decoupling statuses is between 80% and 97%.

Four undesirable decoupling status will emerge if some countries experience negative HDI growth. The worst status is recession decoupling. A significant decline in HDI has accompanied the reduction in carbon emissions per capita. There were more countries with recessive decoupling in the first period. The main countries are the Russian Federation, Cuba, and Mongolia, which were more affected by the collapse of the Soviet Union. Other countries are in Sub-Saharan Africa. In the later periods, the countries with recessive decoupling were mainly in Sub-Saharan Africa, the Middle East, and Latin America. Based on the experience of most countries, the decoupling status between carbon emissions per capita and HDI experienced an evolutionary trend: expansive negative decoupling, expansive coupling, weak decoupling, and strong decoupling. Strong decoupling in the early period is unsuitable for countries with extremely low HDI, which requires vast investments later. However, there has been a shift from recessive decoupling to expansive coupling and expansive negative decoupling for some countries with extremely low HDI (Rwanda, Tanzania, and Zambia).

Overall, the proportion of countries with strong decoupling and weak decoupling continues to rise, while the proportion of expansive coupling and expansive negative decoupling is declining. Although the relationship between HDI and energy consumption has entered a new phase in most countries, the status of strong and weak decoupling is unstable. Strong decoupling and weak decoupling carry the risk of shifting to other status of decoupling. In addition, economic, social, political, and resource factors vary among countries, and the unbalanced trend of economic growth and carbon emissions will

continue in the future. Solving these problems requires the concerted efforts of the whole world.

**Figure 3.** Spatial distribution of the decoupling index over six periods. Source: Authors' own construction using Excel. (**a**) 1990–1995; (**b**) 1995–2000; (**c**) 2000–2005; (**d**) 2005–2010; (**e**) 2010–2015; (**f**) 2015–2019.

Only three countries (the Czech Republic, Germany, and the United Kingdom) have sustained positive decoupling over the six periods. They reduce carbon emissions while maintaining relatively high HDI. Five countries (Mauritius, Panama, South Korea, Turkey, and Uruguay) with high HDI are expansive negative decoupling over at least five periods. Their growth rate of carbon emissions per capita is larger than their HDI. These countries emit more carbon dioxide without a corresponding increase in the HDI.

#### *4.3. Evolution and Two-Dimensional Analysis of Carbon Emission Performance*

#### 4.3.1. Results of Carbon Emission Performance Measurement from 1990–2019

Figure 4 shows the evolution of countries' carbon emission performance in 1990, 1996, 2002, 2008, 2014, and 2019. As shown in Figure 4, there are more and more countries with dark colors. The carbon emission performance of countries is gradually increasing.

**Figure 4.** Evolution of carbon emission performance, 1990–2019. Source: Authors' own construction using Excel. (**a**) 1990; (**b**) 1996; (**c**) 2002; (**d**) 2008; (**e**) 2014; (**f**) 2019.

Countries showing significant increases in carbon emissions performance (50% above from 1990 to 2019) include the Democratic Republic of Congo, Mozambique, Niger, Rwanda, and Sierra Leone, all in sub-Saharan Africa. The most considerable increase was in Rwanda, where CEP increased by 99.89%. In 1990, the total and carbon emissions per capita in Rwanda were 0.52 million tons and 0.0719 tons, respectively. The total and per capita emissions in 2019 were 1.19 million tons and 0.094 tons in Rwanda, respectively. The carbon emission index increased by only 9.69%, whereas the HDI doubled by 119.26%. This result is consistent with the transition of Rwanda's decoupling status. Rwanda started with a deficient level of economic and human development. In the long term, economic growth has accelerated, life expectancy has improved significantly, and primary education has developed significantly in Rwanda. Thus, the increase in HDI results from the overall improvement in all aspects.

Carbon emissions performance declined in 14 countries. The most significant decrease in CEP was recorded in Laos. Total carbon emissions were 0.51 million tons, and carbon emissions per capita were 0.1205 tons in 1990. Total and carbon emissions per capita were 33.93 million tons and 4.7328 tons in 2019, respectively. The carbon emission index of Laos increased by 72.15%, and the human development index increased by 51.35%. The growth of HDI is less than that of carbon emissions. Laos achieved an increase in HDI at the cost of environmental pollution.

The carbon emission performance for the six years was a normal distribution. The majority of countries (73–96%) are at the middle intervals (0.85–1.449). Only a tiny proportion of countries (4–27%) are at the two ends of the interval. The number of countries with CEP at 0.85–1.149 is the greatest, but the proportion has gradually decreased, from 63.64% in 1990 to 53.44% in 2019. The proportion of countries with CEP at 1.15–1.449 gradually increased, rising from 9.79% in 1990 to 42.33% in 2019. The proportion of countries with CEP at 0.55–0.849, mainly in South Asia and Sub-Saharan Africa, fell from 25.17% in 1990 to 2.12% in 2019. Carbon emission performance has been significantly improved, and there are fewer and fewer countries with low carbon emission performance. Although the proportion of countries with CEP at 1.45–1.749 increased to 2.12% in 2019, the number of countries with high CEP is still minimal. There are only four countries (Kiribati, Liechtenstein, Micronesia, and Tonga) with CEP at 1.45–1.79, and their CEP is relatively stable.

#### 4.3.2. Two-Dimensional Analysis Based on HDI and Carbon Emission Performance

The HDI ranges from 0 to 1. The UNDP considers countries with HDI above 0.8 to be very high in human development. The country with the highest CEP is Liechtenstein, with CEP at 1.7471 in 2019. Referring to the analysis method of Zhu and Liu, the upper limit of CEP is set as 2 [18]. According to the classification of HDI, countries with CEP at the top 20% (above 1.6) are regarded as countries with high CEP. The two-dimensional analysis of HDI and CEP in 2019 was performed. Countries can be divided into four types based on the set thresholds: (a) high HDI and high CEP; (b) high HDI and low CEP; (c) low HDI and low CEP; (d) low HDI and high CEP.

As shown in Figure 5, only Liechtenstein has a high HDI and high CEP, which is an ideal state. There are 65 countries with high HDI and low CEP, such as the United States, Norway, Canada, and other developed countries, mainly centered in Europe, Central Asia, and North America. These countries with a very high HDI enjoy better economic and social development. Development should focus on maintaining economic growth while reducing environmental pollution. These countries need to improve emission reduction technologies, develop new renewable energy sources, and optimize industrial structures to achieve energy conservation and emission reduction. The number of countries with low HDI and low CEP is the largest, over 65%, and includes developing countries such as China and India, and least-developed countries, such as Lesotho and Yemen. These countries are in the stage of comprehensive development. They not only need to meet the energy needs to promote the economy and society, but also need to control carbon emissions, and face the dual challenges of economic and social development and environmental sustainability. There are currently no countries with low HDI and high CEP. The fourth type will not appear based on the exponential relationship between carbon emissions per capita and HDI.

#### *4.4. Decomposition of Drivers*

As known from Section 4.2, only the Czech Republic, Germany, and the United Kingdom have maintained continuous strong decoupling over the six periods. The status of decoupling is unstable for most countries, and it is a challenge to sustain strong decoupling for a long time. Therefore, it is necessary to understand the drivers of the decoupling index and CEP of exemplary countries. In addition to the three countries with continuous positive decoupling, the carbon emission paths of five countries with expansive negative decoupling were studied. However, the data for Mauritius, Panama, and Uruguay were not available, and so the decomposition analysis is only for South Korea and Turkey, to explore why they have not yet achieved strong decoupling.

**Figure 5.** Two-dimensional analysis based on HDI and CEP in 2019. Source: Authors' own construction using Excel.

#### 4.4.1. Decomposition of the Decoupling Index

Combining the Tapio decoupling model with the LMDI decomposition model, the decoupling index between carbon emissions per capita and HDI can be understood from the carbon intensity of fossil energy effect, energy structure effect, energy intensity effect, and economic development effect. The smaller decoupling index indicates the lower dependence of HDI growth on energy consumption for strong decoupling status. The strong decoupling index gets smaller and smaller only when the four effects are negative. That is, the contribution rate of the four effects is positive. As shown in Table 2, the decoupling index is becoming smaller, and the strong decoupling effect is increasing in the Czech Republic, Germany, and the UK. The main factors influencing the decoupling index are the same for all three countries. The main contributing factor is the energy intensity effect, indicating high energy use efficiency. On the other hand, the economic development effect is the main inhibiting factor, indicating an increase in carbon emissions in daily life with the popularization of energy-consuming consumer goods. But the contribution of the economic development effect tends to decrease.

In order to transform from expansive negative decoupling to strong decoupling, South Korea and Turkey need to make the decoupling index as small as possible under the premise of positive HDI growth. In other words, all four driving factors are guaranteed to be negative, i.e., the contribution rate is positive. As seen in Table 1, the decoupling index for Korea has decreased in a fluctuating way but has not yet reached weak and strong decoupling. For Korea, the carbon intensity of fossil energy was the main factor contributing to the decline in the decoupling index during the 1990–1995 and 1995–2000 periods. The main factors in the last four periods are transformed into energy intensity. Economic development effects have been the primary inhibitor to transition to weak and strong decoupling for South Korea. Similarly, Turkey's decoupling index has declined, even to the threshold of strong decoupling in the 2015–2019 period. However, the main factors affecting the decline in the decoupling index are different in each period, indicating that some efforts have been made to optimize energy consumption structure and improve energy efficiency. The factor inhibiting the decrease of the decoupling index is always the effect of economic development.

**Table 2.** Decomposition of the decoupling index for the example countries.

Source: Authors' own calculation.

#### 4.4.2. Decomposition of the CEP

As shown in Table 3, the carbon productivity effect, economic development effect, and population effect are positive in most periods for these five countries, indicating that the three effects positively impact the improvement of CEP. But the role of the population effect is relatively small. In most periods, the welfare effect, the carbon intensity of fossil fuels effect, the energy structure effect, and the energy intensity effect are negative. The welfare effect is negative, which means the low CEP is generated by unit economic output growth. The latter three effects are related to the energy utilization, energy structure, and efficiency of fossil fuels, and the negative effects indicate the need to adjust the energy structure further and reduce energy consumption. However, the main influencing factors are different for different countries.

The CEP of the Czech Republic, Germany, and the United Kingdom are primarily similar and gradually increased from around 0.95 in 1990 to 1.2 in 2019, with a growth rate of around 25%. The main influencing factors for the three countries are also the same. The carbon productivity effect is the primary driver of the improvement in CEP. The higher the carbon productivity, the greater the economic output per unit of carbon emissions, representing the progress of emission-reduction technologies in economic development. The main driving factor that inhibits CEP is the energy intensity effect, which indicates that economic growth is highly dependent on primary energy, and new energy needs to be developed to improve energy efficiency.

**Table 3.** Decomposition of CEP for the example countries.

Source: Authors' own calculation.

The CEP of South Korea was 0.95 in 1990. It was 1.09 in 2019, an increase of only 14.74%. The CEP of Turkey was 0.8 in 1990, a low starting point, and then increased significantly to 1.05, at an increase of 31.25%. The main drivers in the two countries are slightly different. In South Korea, it is clear that the most significant positive effect of increasing CEP is the economic development effect, and the most significant negative effect is the welfare effect. Other effects are not significant. The economic development and the living standard in South Korea have greatly improved, which has played a leading role in the improvement of CEP. However, South Korea's economic growth has not promoted the overall welfare improvement very well, resulting in a potent inhibitory effect on the welfare effect. Although the main positive and negative effects in Turkey are also economic development effects and welfare effects, respectively, the contribution rate is not significantly prominent. The facilitation of carbon productivity effects and the inhibition of energy intensity effects also play a role. As a result, CEP in Turkey has been greatly improved. Overall, South Korea and Turkey should focus on energy efficiency and improve the development of emission-reduction technologies.

#### **5. Conclusions and Recommendations**

#### *5.1. Conclusions*

This paper examines the decoupling between carbon emissions and HDI and the CEP. Combining the Tapio decoupling model and the LMDI decomposition model, it is found that the decoupling status and CEP of various human development countries are different. Due to data availability, only the drivers of the decoupling index and CEP for the example countries are decomposed. The main findings of the study are as follows:

#### *5.2. Recommendations*

Some targeted recommendations need to be put forward to achieve decoupling between carbon emissions and HDI and improving CEP.

**Author Contributions:** D.L. conceptualized and designed the research. T.S. drafted the paper and conducted the data interpretation. X.W. and J.L. collected data and obtained preliminary data results. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper is supported by the Key Program of the National Philosophy and Social Science Foundation of China (Grant No. 21ATJ003), Innovation Team of Philosophy and Social Sciences in Henan Colleges and Universities (2017–CXTD-07), Major Projects in Basic Research of Philosophy and Social Sciences in Henan Colleges and Universities (2019-JCZD-002), and Humanities and Social Sciences Project of Henan Educational Administration (2022-ZZJH-416).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The raw data presented in this study are openly available on some platforms of international organizations, reference numbers 1, 26, 27, 28 and 29.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

8. Tapio, P. Towards a theory of decoupling: Degrees of decoupling in the EU and the case of road traffic in Finland between 1970 and 2001. *Transp. Policy* **2005**, *12*, 137–151. [CrossRef]

### *Article* **Study on the Associations between Meteorological Factors and the Incidence of Pulmonary Tuberculosis in Xinjiang, China**

**Chunjie Gao 1,†, Yingdan Wang 1,†, Zengyun Hu 2,3, Haiyan Jiao 1 and Lei Wang 4,\***

**Abstract:** Pulmonary tuberculosis (PTB) has been a major threat to global public health. The association between meteorological factors and the incidence of PTB has been widely investigated by the generalized additive model, auto-regressive integrated moving average model and the distributed lag model, etc. However, these models could not address a non-linear or lag correlation between them. In this paper, a penalized distributed lag non-linear model, as a generalized and improved one, was applied to explore the influence of meteorological factors (such as air temperature, relative humidity and wind speed) on the PTB incidence in Xinjiang from 2004 to 2019. Moreover, we firstly use a comprehensive index (apparent temperature, AT) to access the impact of multiple meteorological factors on the incidence of PTB. It was found that the relationships between air temperature, relative humidity, wind speed, AT and PTB incidence were nonlinear (showed "wave-type ", "invested Utype", "U-type" and "wave-type", respectively). When air temperature at the lowest value (−16.1 ◦C) could increase the risk of PTB incidence with the highest relative risk (RR = 1.63, 95% CI: 1.21–2.20). An assessment of relative humidity demonstrated an increased risk of PTB incidence between 44.5% and 71.8% with the largest relative risk (RR = 1.49, 95% CI: 1.32–1.67) occurring at 59.2%. Both high and low wind speeds increased the risk of PTB incidence, especially at the lowest wind speed 1.4 m/s (RR = 2.20, 95% CI: 1.95–2.51). In particular, the lag effects of low and high AT on PTB incidence were nonlinear. The lag effects of extreme cold AT (−18.5 ◦C, 1st percentile) on PTB incidence reached a relative risk peak (RR = 2.18, 95% CI: 2.06–2.31) at lag 1 month. Overall, it was indicated that the environment with low air temperature, suitable relative humidity and wind speed is more conducive to the transmission of PTB, and low AT is associated significantly with increased risk of PTB in Xinjiang.

**Keywords:** pulmonary tuberculosis; penalized distributed lag non-linear model; meteorological factors; apparent temperature; cumulative risk

### **1. Introduction**

As one of the widely distributing and potentially fatal infectious diseases, tuberculosis remains one of the top ten causes of mortality worldwide. According to the Global Tuberculosis Report of the World Health Organization in 2020 [1], there were approximately 9.96 million cases and 1.42 million tuberculosis-related deaths worldwide in 2019, of which about 5.94 million were new PTB cases. China accounted for 8.6% of the world's new tuberculosis cases in 2019, ranking third among 30 countries with a high burden of tuberculosis. The incidence of PTB in Xinjiang is relatively serious, which has always been one of the highest in China, and the reported incidence of PTB in 2019 is about three times the national level.

**Citation:** Gao, C.; Wang, Y.; Hu, Z.; Jiao, H.; Wang, L. Study on the Associations between Meteorological Factors and the Incidence of Pulmonary Tuberculosis in Xinjiang, China. *Atmosphere* **2022**, *13*, 533. https://doi.org/10.3390/atmos 13040533

Academic Editor: Kai-Jen Chuang

Received: 24 February 2022 Accepted: 24 March 2022 Published: 28 March 2022

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

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

PTB is a chronic respiratory disease caused by *Mycobacterium tuberculosis* (*M. tuberculosis*). A large number of lipids in *M.tb* resist the multiplication of lysozymes and damage macrophages. Hence, it can inhibit the intracellular bactericidal mechanism and cause an inflammatory response in the lungs [2]. Susceptible people can be infected through breathing in dust with *M.tb* and inhaling the droplet nucleus released into the air by actively infectious PTB coughing or sneezing [3]. It is revealed that meteorological factors (such as air temperature, relative humidity and wind speed) may indirectly affect *M.tb* transmission in the environment [4–7]. For instance, within a certain range of air temperature and relative humidity, droplets containing *M.tb* are more likely to be evaporated in the air to form into certain diameters that can be suspended in the air for a longer time (especially, droplets containing *M.tb* can grow and reproduce at 35~37 ◦C, and it can survive for 4~5 years in −8~−6 ◦C), to be easily inhaled into the body by susceptible persons [5]. Moreover, changes in air temperature can also affect the human body's physiological response to toxic agents and retard the clearance rate of *M.tb* [6], and activate or inhibit the development of PTB by influencing blood pressure [7]. Suitable humidity can promote the growth and reproduction of *M.tb* and increase the time for *M.tb* to float in the air [8]. When *M.tb* adheres to dust, it is easily affected by the wind speed, i.e., the higher the speed of the wind, the wider spread of *M.tb* [9].

The annual PTB incidence in Xinjiang is still higher, which may be closely related to its unique geographical location and climatic conditions [10]. Xinjiang, located at the northwest border of China, has a temperate continental climate, with long winters, sparse precipitation, dry climate, frequent wind and sand activities, and frequent snowfall [11]. With an annual average temperature of 11 ◦C, Xinjiang is one of the regions with a lower annual average temperature in China. The heating period in Xinjiang is from October to March of the next year, and the air pollution is serious in this period. *M.tb* is so small that normal air currents can keep the pollution particulates containing *M.tb* airborne and transport them through rooms or some buildings, which may increase the risk of PTB development [3,5,12]. In addition, the snowfall in winter increases the humidity and makes the droplet nucleus ejected by actively infectious PTB stay in the air longer [13], which may create favorable conditions for the transmission of PTB. Therefore, it is of great significance to study the relationship between meteorological factors and the PTB incidence in Xinjiang.

There have been various statistical models to characterize the relationship between meteorological factors and PTB incidence. For instance, a generalized additive model [14] was used to quantitatively evaluate the effects of meteorological factors on the risk of pulmonary tuberculosis in Jiangsu Province. The results illustrated that an environment with low temperature, relatively high wind speed, and low relative humidity is conducive to the transmission of PTB. Li et al. [15] employed an auto-regressive integrated moving average model to exhibit the best predicting performance of PTB incidence by incorporating meteorological factors. However, some results [8,9,15,16] demonstrated there is a non-linear and lag correlation between meteorological factors and the incidence of PTB. The generalized additive model somewhat ignores the collinearity among different lag days [4,17], and the auto-regressive integrated moving average model fails to solve the collinearity problem and ignores the lag effects [12,18]. Fortunately, the distributed lag non-linear model (DLNM) proposed by Armstrong [19] can effectively evaluate the nonlinear and lag relationship between them. For example, Wu et al. [20] investigated the cold and hot effects on mortality at different lags in four subtropical cities. Yang et al. [21] evaluated the association between meteorological factors and the mumps incidence in Guangzhou. Subsequently, a penalized framework DLNM (P-DLNM) was put forward in [22], which introduces a penalized framework to the exposure and lag dimensions in the DLNM, to overcome the complex fitting and the lack of general standards in the selection of basic functions, the number of nodes, the maximum number of lag days and the optimal model in the DLNM. It can control the smoothness of the basis function, and effectively reduce the total number of degrees of freedom to make itself more robust.

Another concerning issue is how the combination of some meteorological factors influenced the incidence of PTB. AT, firstly proposed by Stead-man [23], is a comprehensive index to measure the interaction of air temperature, wind speed and relative humidity, which more objectively represents the actual perception of air temperature. Many research works [15,16,24] evaluated the relationship between meteorological factors and human health by applying AT, such as acute coronary syndromes [15], stroke [16] and acute excessive drinking [24]. However, there are few studies that quantified the impact of meteorological factors and AT on PTB incidence in Xinjiang. Therefore, in this paper, a penalized framework DLNM model is used to analyze the impact of meteorological factors (such as air temperature, wind speed, and relative humidity) on PTB incidence, and AT, is applied to quantify the comprehensive impact of multiple meteorological factors on PTB incidence in Xinjiang from 2004 to 2019, which provides a theoretical basis for the prevention and control of PTB in Xinjiang.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Xinjiang, located in Northwest China (73◦40 –96◦18 E, 34◦25 –48◦10 N), has a special topography named as two basins lie in between three mountains, with an average altitude of about 1000 km. From 2004 to 2019, the average annual temperature in Xinjiang was 10.1 ◦C, the average annual precipitation was 110.0 mm [11], the average annual wind speed was 2.3 m/s, the average annual air pressure was 898 hPa and the average annual relative humidity was 48.3% [25].

#### *2.2. Meteorological and PTB Data*

The monthly PTB cases in Xinjiang from 2004 to 2019 were obtained from Public Health Scientific Data Sharing Center (http://www.phsciencedata.cn/, accessed 30 December 2021)and Health Commission of Xinjiang Uygur Autonomous Region (http://wjw.xinjiang. gov.cn/, accessed 8 January 2022).

The monthly average values of meteorological indicators from 2004 to 2019 were derived from 54 stations in Xinjiang (Supplementary Figure S1) from the China Meteorological Data Sharing Center (http://data.cma.cn, accessed 21 November 2021), including average air temperature (◦C), average precipitation (mm), average wind speed (m/s), average air pressure (hPa), average relative humidity (%), average sunshine duration (hours/month). AT, is calculated as follows [24]:

$$AT = T + 0.33 \times e - 0.70 \times W - 4.00,\tag{1}$$

$$
\omega = \frac{RH}{100} \times 6.105 \times \exp\left(17.27 \times \frac{T}{237.7 + T}\right),
\tag{2}
$$

where *T* represents air temperature (◦C); *e* represents water vapor pressure (hPa); *RH* is relative humidity (%); and *W* represents wind speed (m/s).

In this paper, based on percentile range the AT was divided into seven categories for analysis of the associations between low AT, high AT and PTB incidence: extreme cold (≤1st percentile), cold (1st–5th percentiles), mild cold (5th–25th percentiles), comfortable (25th–75th percentiles), mild heat (75th–95th percentiles), heat (95th–99th percentiles) and extreme heat (≥99th percentile), respectively. In addition, the influences of low AT (including extreme cold, cold and mild cold) and high AT (including extreme heat, heat, and mild heat) on PTB incidence were investigated [24].

#### *2.3. Spearman's Rank Correlation*

In order to measure the nonlinear relationship between meteorological variables and PTB incidence, Spearman's rank correlation analysis was used [26]. The correlation coefficient *rs* denotes the strength of an association between two variables, which can be calculated as follows [27]:

$$r\_s = 1 - \frac{6\sum d\_i^2}{n(n^2 - 1)}$$

where *n* is the count of variables and *d* is the rank difference between each pair of variables. The heatmap can be used to represent the values of the correlation coefficient between multiple variables. In the heatmap, the darker the color grid, the greater the Spearman's rank correlation coefficient, and the blank space is usually used to represent two variables without significant correlation.

#### *2.4. Statistical Method*

DLNM, proposed by Gasparrini and Armstrong in 2010 [19], is widely used to evaluate the relationship between meteorological factors and health effects, which is a nonlinear model reflecting the exposure-lag-response relationship by selecting basis functions for exposure-response and lag-response and to form a cross-basis function [28].

Define, *f*(*xt*; *β*) = *Zt Tβ* as the exposure-response function, where *Zt* is a matrix with dimensions *n* × *νx* (*n* represents the length of the time series, *νx* represents the degree of freedom), and *β* is the shape parameter of *Zt T*. Define, *qx*,*t* = *xt*−-0 ,..., *xt*−-,..., *xt*−*L T* is lagged effects of the exposure *xt*, where -0 and *L* are the minimum and maximum lags, respectively, i.e., the lag structure can be formulated as [28] - = [-0,..., *L*] *T*. The cross-basis function *s*(*xt*, η) is constructed by:

$$s(\mathbf{x}\_{t\prime}\boldsymbol{\eta}) = \left(\mathbf{1}\_{L-\ell\_0+1}^T \boldsymbol{A}\_{\mathbf{x},t}\right) \boldsymbol{\eta} = \boldsymbol{w}\_{\mathbf{x},t}^T \boldsymbol{\eta}\_{\prime\prime}$$

where *wT x*,*t* is a cross-basis matrix transformed from *Ax*,*t* (with dimensions *νx* × *ν*-), and *η* is the coefficient of *wT x*,*t*. *Ax*,*t* is computed by a row-wise Kronecker product between the two basis matrices *R* (with dimensions (*L* − -0 + 1) × *νx*) and *C* (with dimensions (*L* − -0 + 1) × *ν*-) that obtained by applying basis transformations to vectors *qx*,*t* and -, thus *Ax*,*t* can be expressed as:

$$\mathbf{A}\_{\mathbf{x},\mathbf{t}} = \left(\mathbf{1}\_{\nu\_{\ell}}^{T} \otimes \mathbf{R}\_{\mathbf{x},\mathbf{t}}\right) \odot \left(\mathbf{C} \otimes \mathbf{1}\_{\nu\_{\nu\_{\mathbf{x}}}}^{T}\right).$$

where *1T ν* is a unit vector with length *ν*-, ⊗ and  denote the Kronecker and Hadamard products, respectively [29].

The distributed lag non-linear model with a penalized framework (P-DLNM) applies varying degrees of penalties in exposure and lag dimensions to smooth the surface of exposure-lag-response, the log-likelihood function *lp*(*η*, *γ*, *λ*) of P-DLNM is:

$$l\_p(\boldsymbol{\eta}, \boldsymbol{\gamma}, \boldsymbol{\lambda}) = l(\boldsymbol{\eta}, \boldsymbol{\gamma}) - \frac{1}{2} \boldsymbol{\eta}^T \left(\lambda\_x \left(\mathbf{S}\_x \otimes \mathbf{1}\_{\boldsymbol{\nu}\_\ell}^T\right) + \lambda\_\ell \left(\mathbf{1}\_{\boldsymbol{\nu}\_\ell}^T \otimes \mathbf{S}\_\ell\right)\right) \boldsymbol{\eta}\_\ell$$

where the penalty matrices *Sx* and *S* are separately obtained by introducing the penalties term to the base matrices *R* and *C*, *λx* and *λ* represent the penalty parameters of *Sx* and *S*-, respectively. *η* is the coefficient for the cross-basis, restricted by *Sx*, *S* and the penalty parameter *λ* = [*λx*, *λl*] *T*. In this paper, the ridge penalty is used to add the penalty term for the lag dimension and control the coefficient of *S* to shrink towards the zero value at different lag times [29]. The penalties of lag-response dimension are based on cubic regression splines, i.e., CR smoother, which can be expressed as *S*- = *Pν*- , here *Pν* is a pre-specified diagonal matrix with the weight *p*.

Based on the generalized additive model and combined the cross-basis functions with penalties, meanwhile, controlled the long-term trend of time and other confounding factors, the association between meteorological factors and PTB incidence was investigated by the following model in this paper:

$$\ln(\lambda\_I) = a + \sum\_{i=1}^{4} s\_i(\mathbf{x}\_i, \eta\_i) + ns(MF, df) + ns(time, df) + \gamma S + \varepsilon AR + \theta SF + \tau PC\_{\tau}$$

where *λt* is the case number of PTB on month *tth*; *α* is the intercept; are the cross-basis functions with coefficients *ηi* (*i* = 1, . . . , 4), *xi* are meteorological variables (such as air temperature, wind speed, relative humidity and AT); *ns*(·) is a natural cubic spline function, which can be used to control the long-term trends (time with the degree of freedom 10, *d f* = 10) and other meteorological factors (air pressure, precipitation and sunshine duration with the degrees of freedom 3, 3 and 1, respectively). *S* is a categorical variable with coefficient *γ* to control for the seasonality (spring: Mar.-May., summer: Jun.-Aug., autumn: Sep.-Nov. and winter: Dec.-Feb.); *AR* is the autoregressive term with coefficient *ε* to correct for autocorrelation in the residuals [30]; is a binary variable with coefficient *θ* to control for the potential influences on delayed medical treatment of PTB patients caused by the Spring Festival [12]; *PC* is a binary variable with coefficient *τ* to avoid the impact of PTB transmission caused by immigrants entering Xinjiang to pick cotton from August to October every year.

In this paper, taking the median of air temperature, wind speed, relative humidity and AT as references, four cross-basis functions (defined by an exposure-response relationship with four equally spaced intervals knots and a lag-response relationship with four equal space log-knots) were built to more flexibly reflect the 3-D relationship between meteorological factors and the incidence of PTB in Xinjiang.

#### *2.5. Parameter Estimation*

To quantify the risk between meteorological factors and PTB incidence, the relative risks (*RR*) of specific exposure and lag time were used [22]. The values of *RR* and standard error (*SE*) of specific exposure level *xp* are as follows:

$$RR\_{x\_p} = \exp\left(\hat{\mathcal{J}}\_{\mathbf{x}\_p}\right) = \exp\left(A\_{\mathbf{x}\_p}\hat{\eta}\right)$$

$$SE\_{x\_p} = \exp\left(\sqrt{V\left(\hat{\mathcal{J}}\_{\mathbf{x}\_p}\right)}\right) = \exp\left(\sqrt{A\_{\mathbf{x}\_p}V(\hat{\eta})A\_{\mathbf{x}\_p}^T}\right)$$

where *Axp* denotes the cross-basis matrix with *x* = *xp*; *V*(·) denotes the variance of a random variable. Similarly, the lag-specific risk *RRp* and *SEp* can be estimated. In addition, the overall cumulative effect (*RRc*) of such exposure history, with associated *SEc* could be computed with:

$$RR\_{\mathfrak{c}} = \exp\left(\hat{\mathfrak{P}}\_{\mathfrak{c}}\right) = \exp(\mathbf{W}^{p}\hat{\mathfrak{p}}),$$

$$SE\_{\mathfrak{c}} = \exp\left(\sqrt{V\left(\hat{\mathfrak{P}}\_{\mathfrak{c}}\right)}\right) = \exp\left(\sqrt{\mathbf{W}^{p}V(\hat{\mathfrak{p}})\mathbf{W}^{pT}}\right).$$

where *Wp* denotes the cross-basis matrix.

#### *2.6. Sensitivity Analysis*

To assess the robustness of the model, a sensitivity analysis was conducted to examine the influence of *d f* of parameters and maximum lag time on the overall cumulative effect of AT. In this paper, *d f* of time in the model was varied from 10, 11 to 12, respectively. The *d f* of air pressure and precipitation in the model were changed between 3 and 5, respectively. The *d f* of sunshine duration was altered between 1 and 3, respectively. The maximum lag time of AT was also set to 8, 10 and 12, respectively.

All statistical analyses were performed by carrying out R software version 4.0.5.

#### **3. Results**

#### *3.1. Descriptive Statistics of PTB Cases and Meteorological Factors*

During the study period from 2004 to 2019 (192 months), a total of 668,753 PTB cases were included, with an average of 3427 cases per year and a maximum of 74,549 cases in 2018. The case number of PTB in Xinjiang had an obvious seasonal pattern and showed an increasing trend, with a peak from January to April (see Figure 1 and Supplementary Figure S2).

**Figure 1.** Time series of monthly PTB cases and meteorological factors in Xinjiang from 2004 to 2019. Abbreviations: Temp, air temperature (◦C); AT, apparent temperature (◦C); SD, sunshine duration (h); RH, relative humidity (%); AP, air pressure (hPa); Pre, precipitation (mm); WS, wind speed (m/s). Note: The blue dotted line represents the median.

There were seasonal fluctuations and periodic trends of meteorological factors in Xinjiang, roughly showing the variation of single peak and single valley (see Figure 1 and Supplementary Figure S3). Air temperature, AT and sunshine duration had similar seasonal patterns, with higher values occurring from April to October (Seasonal index >1 indicates obvious seasonal variation). The peak of relative humidity and air pressure occurred from October to February of the next year, whereas the trough appeared from March to September. A seasonal peak of precipitation emerged from May to September and a trough occurred from October to March of the next year. Wind speed had a seasonal trend with a peak from March to August and a valley from September to February of the next year. In addition, the median of AT, air temperature, relative humidity, air pressure, precipitation, wind speed and sunshine duration were 7.8 ◦C, 11.3 ◦C, 45.6%, 898.2 hPa, 8.0 mm, 2.3 m/s and 251.0 h, respectively (see Figure 1).

The descriptive statistics of meteorological factors and the case number of PTB were shown in Table 1. The monthly average of AT was 5.7 ◦C, which was 3.5 ◦C lower than the air temperature. The distributions of AT in the1st, 5th, 25th, 75th, 95th and 99th percentiles were −18.5 ◦C, −15.1 ◦C, −6.3 ◦C, 7.8 ◦C, 18.0 ◦C, 23.3 ◦C and 24.3 ◦C, respectively.

**Table 1.** Descriptive statistics for monthly PTB cases and meteorological factors in Xinjiang from 2004 to 2019.

Abbreviations: Temp, air temperature (◦C); AT, apparent temperature (◦C); SD, sunshine duration (h); RH, relative humidity (%); AP, air pressure (hPa); Pre, precipitation (mm); WS, wind speed (m/s).

#### *3.2. Spearman's Rank Correlation Analysis*

The Spearman's rank correlation results between monthly PTB cases and meteorological factors in Xinjiang from 2004 to 2019 were provided in Figure 2. AT (*rs* = −0.25), air temperature (*rs* = −0.23), sunshine duration (*rs* = −0.17) and precipitation (*rs* = −0.14) were negatively correlated with the incidence of PTB, and wind speed (*rs* = 0.15) is positively correlated with it. There was no significant correlation between relative humidity, air pressure and the incidence of PTB. There were high correlations among meteorological factors, especially between air pressure and AT (*rs* = −0.94). The *p*-value < 0.05 was considered statistically significant.

**Figure 2.** Spearman's rank correlation results between the cases number of PTB and meteorological factors in Xinjiang from 2004 to 2019. Abbreviations: Temp, air temperature (◦C); AT, apparent temperature (◦C); SD, sunshine duration (h); RH, relative humidity (%); AP, air pressure (hPa); Pre, precipitation (mm); WS, wind speed (m/s); *rs*, Spearman's rank correlation coefficient.

#### *3.3. The Influences of Air Temperature on the Incidence of PTB*

Taking the overall median of monthly average air temperature (11.3 ◦C) as a reference, it was obtained that the correlation between average air temperature and PTB incidence was nonlinear at 12 lag months for the general population, and the largest exposed and delayed cumulative effects occurred at −16.1 ◦C and lag 0 months, respectively (see Figure 3A). The overall exposure-response curve between air temperature and the incidence of PTB took

on a U-type when the air temperature below 11.3 ◦C, with a maximum value, occurred at −16.1 ◦C (RR = 1.63, 95% CI: 1.21–2.20), whereas an inverted U-type when the air temperature was higher than 11.3 ◦C with a minimum value appeared at 26 ◦C (RR = 0.12, 95% CI: 0.06–0.20), as shown in Figure 3B. The association between extreme air temperature and the risk of PTB at specific lag months was shown in Figure 3B and Table S1. The lag 0–1 and lag 11–12 months could significantly increase the risk of PTB at extremely low air temperature −11.9 ◦C (2.5th percentiles), with a peak at lag 0 (RR = 1.59, 95% CI: 1.53–1.66). The cumulative risk at lag 12 months was largest (RR = 1.19, 95% CI: 1.14–1.24) at the extremely high air temperature 25 ◦C (97.5th percentiles). The relationship of air temperature-lag-PTB was evaluated by using the heatmap, as shown in Figure 3C.

**Figure 3.** The exposure-lag-response correlations between meteorological factors and PTB incidence in Xinjiang from 2004 to 2019. (**A**–**C**) represent the 3–D graphs of the relationships between air temperature, relative humidity, wind speed and PTB incidence, respectively. (**D**–**F**) represent the overall cumulative relative risks of air temperature, relative humidity, wind speed on PTB incidence acrosslag 0–12 months, respectively. (**G**–**I**) represent the heatmaps of air temperature, relative humidity, wind speed on PTB incidence, respectively.

#### *3.4. The Influences of Relative Humidity on the Incidence of PTB*

Setting the overall median of monthly average relative humidity (45.6%) as a reference, it was found that the association between relative humidity and PTB incidence was nonlinear at 12 lag months for the general population, and the largest exposed and delayed cumulative effects occurred at the relative humidity 59.2% and lag 0 months, respectively (see Figure 3D). The overall exposure-response curve between relative humidity and the PTB incidence appeared an inverted U-type, with a peak at a relative humidity of 59.2% (RR = 1.49, 95% CI: 1.32–1.67), as shown in Figure 3E. The relationship between extreme relative humidity and the risk of PTB at specific lag months was revealed in Figure 3E and Table S1. The lag 0, 10, 11 and 12 months could significantly increase the risk of PTB at the extreme low relative humidity of 32.8% (2.5th percentiles), with a top at lag 0 (RR = 1.05, 95% CI: 1.03–1.07). The lag 0–3 months could significantly increase the risk of PTB at the extreme high relative humidity 68.9% (97.5th percentiles), with a top at lag 0 (RR = 1.15, 95% CI: 1.13–1.17). The association of relative humidity-lag-PTB was estimated by using a heatmap, as shown in Figure 3F.

#### *3.5. The Influences of Wind Speed on the Incidence of PTB*

Using the overall median of monthly average wind speed of 2 min (2.3 m/s) as a reference (the wind speed corresponding to the lowest risk of PTB), it was found that the influence of wind speed on PTB incidence was nonlinear at 12 lag months for the general population, and the largest exposed and delayed cumulative effects occurred at 1.4 m/s and lag 2 months, respectively (see Figure 3G). The overall exposure-response curve between wind speed and the incidence of PTB showed a U-type with a trough that emerged at 2.3 m/s and a maximum relation arose at 1.4 m/s (RR = 2.20, 95% CI: 1.95–2.51), as shown in Figure 3H. The correlation between extreme wind speed and the risk of PTB at specific lag months was shown in Figure 3H and Table S1; the lag 0–9 months were separately significantly associated with an increased risk of PTB effect at the extreme low wind speed of 1.6 m/s (2.5th percentiles), with a maximum value at lag 2 (RR = 1.15, 95% CI: 1.13–1.16). The lag 6–12 were separately significantly related to an increased risk of PTB at the extreme high wind speed 3.1 m/s (97.5th percentiles), with a maximum value at lag0 (RR = 1.14, 95% CI: 1.12–1.15). The correlation of wind speed-lag-PTB was also assessed by using the heatmap as shown in Figure 3I.

#### *3.6. The Effect of AT on PTB Incidence*

Considering the overall median of the monthly AT (7.8 ◦C) as a reference, it was discovered that the connection between AT and the incidence of PTB was nonlinear at 12 lag months for the total population, and the largest exposed and delayed cumulative effects occurred at −20.9 ◦C and lag 0 months (see Figure 4A). As shown in Figure 4B, the overall exposureresponse curve between AT and the incidence of PTB presented a *U*-type when the AT below 7.8 ◦C, with a maximum value occurred at −20.9 ◦C (RR = 2.01, 95% CI: 1.50–2.69), whereas an inverted *U*-type when the AT was higher than 7.8 ◦C with a maximum risk appeared at 12.6 ◦C (RR = 1.29, 95% CI: 1.16–1.44). The extreme cold −18.5 ◦C and mild heat 18 ◦C had no significantly risky effects on PTB incidence, and the cold −15.1 ◦C (RR = 0.59, 95% CI: 0.46–0.76), the mild cold −6.3 ◦C (RR = 0.22, 95% CI: 0.18–0.27), the heat 23.3 ◦C (RR = 0.58, 95% CI: 0.38–0.90) and the extreme heat 24.3 ◦C (RR = 0.51, 95% CI: 0.32–0.82) had significantly protective effects on PTB incidence. The connection of AT-lag-PTB was also observed by using the heatmap, as shown in Figure 4C.

It is noteworthy that the lag effects of low AT presented a *V*-shape with the highest risk at lag 0, whereas the lag effects of high AT presented an *N*-shape with the highest risk at lag 12 (see Figure 5). As shown in Figure 5, lag 0 and lag 1 could significantly increase the incidence of PTB at the extreme cold, with the highest risk at lag 0 (RR = 1.81, 95% CI: 1.73–1.90). The cold had a significantly risky effect on increasing PTB incidence at lag 0–1, with the highest risk at lag 0 (RR = 1.60, 95% CI: 1.54–1.67), which was a protective factor for the risk of PTB at lag 2–11. The mild cold had a significantly risky effect on increasing the incidence

of PTB at lag 0 and lag 1, with the highest risk at lag 0 (RR = 1.22, 95% CI:1.19–1.26), which was a protective factor for the risk of PTB at lag 2–12. The mild heat had a significantly risky effect on increasing the incidence of PTB at lag 2 and lag 9–12, with the highest risk at lag 12 (RR = 1.11, 95% CI: 1.09–1.14), which did not significantly influence the incidence of PTB at lag 3. The heat and extreme heat had a significantly risky effect on increasing the incidence of PTB at lag10–12, with the highest risk at lag 12 (RR = 1.18, 95% CI: 1.14–1.22 and RR = 1.20, 95% CI: 1.15–1.25), which did not significantly influence PTB incidence at lag 2 and lag 8–9.

**Figure 5.** The distribution of RR at various lag months with a reference of medium value 7.8 ◦C in Xinjiang from 2004 to 2019.

In addition, the cumulative effects of low and high AT on PTB incidence at different lag months were compared (see Table 2). For the extreme cold of low AT, the cumulative effects were risky from lag 0 to lag 0–8, with a maximum value at lag 0–1 (RR = 2.18, 95% CI: 2.06–2.31). For the cold, the cumulative effects were risky effects from lag 0 to lag 0–4, with a maximum value at lag 0–1 (RR = 1.84, 95% CI: 1.74–1.94). For the mild cold the cumulative

effects were risky from lag 0 to lag 0–2, with a maximum value at lag 0–1 (RR = 1.28, 95% CI: 1.22–1.34). For the high AT, there was no significant risk of the cumulative effects of the mild heat, heat and extreme heat, in particular, the cumulative effects of mild heat had no significant influences from lag 0–10 to lag 0–12.

**Table 2.** The cumulative effects of low and high AT on the incidence of PTB by different lag period in Xinjiang from 2004 to 2019.

#### *3.7. Sensitivity Analysis*

The result of sensitivity analysis showed the model was robust when the *df* were altered for the time trend (*df* = 10–12), air pressure (*df* = 3–5), precipitation (*df* = 3–5) and sunshine duration (*df* = 1–3) (see Figures S4–S6). Changing the maximum lag day into 8, 10 and 12 in the model did not show an obvious difference for fitting the overall effect curve of AT in the model (see Figure S7).

#### **4. Discussion and Conclusions**

In this paper, the effects of meteorological factors (air temperature, wind speed and relative humidity) on PTB incidence in Xinjiang were investigated by using a P-DLNM model, and the relationship between AT and PTB incidence was evaluated. It was found that the correlations between air temperature, relative humidity, wind speed, AT and PTB incidence were nonlinear and lagged.

The overall effect between air temperature and PTB incidence showed a curve fluctuation, and the low air temperature could have a significantly risky effect on PTB incidence, which was consistent with previous studies [4,31,32]. It was reasonable that Xinjiang has a longer winter and thus, people are more susceptible to being infected with *M.tb* because the time of outdoor activities for humans was shorted [31]. Another reason is that the air pollution is serious during the heating period in Xinjiang (from October to March of the next year), then pollution particulates may attach more pathogenic bacteria including *M.tb* [3,12]. Moreover, vitamin D was demonstrated as an important factor that affected people's immune response to resist and remove *M.tb* when someone was during incubation, and with a long winter and low levels of vitamin D, latent infected PTB peoples are sensitive to *M.tb* [32]. The association between high air temperature with PTB incidence is negative, which may be the result that the recombinant strain of *M.tb* may stop growing or even be destroyed when the temperature exceeds 37 ◦C. In addition, sunshine is the main route of vitamin D synthesis. More sunbathing could help humans to enhance immunity so that the risk of transmission of PTB is relatively lower [33]. As a reference of the overall median of monthly average relative humidity (45.6%), it was obtained that the relative humidity was possible to increase the risk of PTB when it was between 44.5% and 71.8%, while it decreased the risk if it ≥ 72%, which is consistent with the results in reference [4]. One explanatory hypothesis for the high relative humidity is a lower chance of infection with

PTB because Xinjiang is an arid/semi-arid region with a dry climate, scarce precipitation and infrequent high relative humidity environment [11]. Another explanatory hypothesis is high relative humidity promotes the production of protective mucus on the surface of the respiratory tract thereby resisting the invasion of *M.tb* [34]. The overall effect of wind speed on the incidence of PTB was a *U*-shape, low and high wind speed are significantly positively correlated with the incidence of PTB, which is different from previous results that high wind speed was positively correlated with PTB [31,35]. It can be seen in the results that high wind speed may affect the transmission of *M.tb*, which may be because the frequent dust weather in Xinjiang makes *M.tb* float in the air for a longer time [9]. Moreover, the lower wind speed may also increase the chance of *M.tb* floating in the air in Xinjiang because the wind speed in winter is lower but the frequency is higher [36].

One of the highlights of this paper is to illustrate that there is a significant association between low AT and PTB incidence in Xinjiang. Generally, people are exposed to multiple meteorological factors in the environment, with low air temperature, wind and high relative humidity the sensations experienced by the human body are usually lower than the air temperature. AT is a comprehensive index combining ambient temperature, wind speed and humidity, which characterizes the physiological experience better than just air temperature alone, then it may be more realistic and objective. In this paper, it was shown that there was a nonlinear and lagged relationship between AT and PTB incidence in Xinjiang, and low AT can significantly increase the risk of PTB. Using the overall median of monthly AT (7.8 ◦C) as a reference, it was discovered that the lag effects of low AT presented a *V*-shape with the highest risk at lag 0, while the lag effects of high AT emerged as an *N*-shape with the highest risk at lag 12. Some specific biological mechanisms can account for why low AT can significantly increase the risk of PTB [37]. In a low air temperature, the human body is more sensitive to cold stimulation than heat stimulation, and the sympathetic nervous system is stimulated if the human body is invaded by cold air, which leads to the contraction of bronchial smooth muscle and excessive response, thus affecting the pulmonary ventilation and the decline of human lung function [38]. Furthermore, the human body's perception of cold and heat is mainly controlled by the thermoregulation system, nervous system and endocrine system [39], *M.tb* is more likely to invade the body in extreme cold weather, thus the human body's autonomous regulation is limited, the defense ability of the immune system is weakened. AT, as a comprehensive index combining air temperature, wind speed and relative humidity can be used to evaluate the impact of meteorological factors on PTB incidence, thus it can better reflect the relationship among human body, environment and disease.

The adjustment of energy structure has improved air quality in Xinjiang, and the incidence of PTB decreased slowly in recent years [25]. However, there is still a long way to go to mitigate the adverse effect of meteorological factors on the incidence of PTB in Xinjiang. Some measures should be implemented to improve the air quality. For example, the government should encourage the use of clean energy during the heating period. Moreover, more trees should be planted to alleviate the harm of dust storms. People should be appealed to wear masks in public places (especially in dusty weather) and enhance their immunity by getting more sunlight (at a favorable temperature) and exercising.

This paper also has some limitations. Firstly, the incidence of PTB is also affected by many variables (such as age, gender and air pollutant concentration) [40]. If these confounding variables are included in the model, the impact of meteorological factors and AT on PTB incidence can be more accurately evaluated. Secondly, other periodic dynamic models also can be applied to investigate the relationship between meteorological factors and PTB incidence.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/atmos13040533/s1, Figure S1: The study area and locations of weather stations; Figure S2. The plot of partial auto-correlation function in the distributed lag non-linear model; Figure S3. Seasonal decomposition of the time series of meteorological factors and the case in-cidence of PTB in Xinjiang from 2004 to 2019. Abbreviations: PTB, pulmonary tuberculosis; Temp, temperature; AT, apparent temperature; SD, sunshine duration; RH, relative humidity; AP, air pressure; Pre, precipitation; WS, wind speed; Figure S4. Sensitivity analysis when altering the degrees of freedom (*df* = 10–12) of time for con-trolling for the long-term trend in the model in Xinjiang from 2004 to 2019; Figure S5. Sensitivity analysis when altering the degrees of freedom (*df* = 3–5) of air pressure and precipitation for controlling for the effect of confounding factors in the model in Xinjiang from 2004 to 2019; Figure S6. Sensitivity analysis when altering the degrees of freedom (1–3) of sunshine duration for controlling for the effect of confounding factors in the model in Xinjiang from 2004 to 2019; Figure S7. Sensitivity analysis when altering the maximum lag periods for 8, 10 and 12 months in the model in Xinjiang from 2004 to 2019; Table S1. Estimated relative risks (95% CI) of pulmonary tuberculosis cases with extremely low temperature (2.5th percentile, −11.9 ◦C) and ex-tremely high temperature (97.5th percentile, 25 ◦C) and extremely low relative humidity (2.5th percentile, 32.8%) and extremely high relative humidity (97.5th percentile, −68.9%) and extremely low wind speed (2.5th percentile, 1.6 m/s) and extremely high wind speed (97.5th percentile, 3.1 m/s) at lagged months in Xinjiang from 2004 to 2019.

**Author Contributions:** Conceptualization, C.G. and L.W.; Data curation, H.J.; Formal analysis, L.W.; Funding acquisition, L.W.; Methodology, C.G. and Y.W.; Software, Y.W. and H.J.; Writing—original draft, C.G., Y.W. and L.W.; Writing—review & editing, C.G., Z.H. and L.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work is supported by the Natural Science Foundation of Xinjiang (Grant No. 2019D01C20), the National Natural Science Foundation of China (Grant Nos. 12061079, 12026236, 12026222), the Shenzhen science and technology innovations Committee (JCYJ20210324101406019) and the Alliance of International Science Organizations (Grant No. ANSO-CR-KP-2021-02).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used and/or analyzed in the current study are publicly available in [the Public Health Scientific Data Sharing Center] at [http://www.phsciencedata.cn/. Accessed 30 December 2021], [the Health Commission of Xinjiang Uygur Autonomous Region] at [http://wjw.xinjiang.gov.cn/. Accessed 8 January 2022] and [the China Meteorological Data Sharing Center] at [http://data.cma.cn. Accessed 21 November 2021].

**Conflicts of Interest:** The authors declare no competing interests.

#### **References**

### *Article* **Impacts of Climate and Environmental Change on Bean Cultivation in China**

**Sidan Li 1,\*, Shibing You 1, Ze Song 1, Li Zhang 1 and Yixuan Liu 2**

**Abstract:** The impact of human-caused environmental pollution and global climate change on the economy and society can no longer be underestimated. Agriculture is the most directly and vulnerably affected sector by climate change. This study used beans, the food crop with the largest supply and demand gap in China, as the research object and established a panel spatial error model consisting of multiple indicators of four factors: climate environment, economic market, human planting behavior and technical development level of 25 provinces in China from 2005 to 2019 to explore the impact of climate environmental changes on the yields of beans. The study shows that: (1) The increase in precipitation has a significant positive effect on bean yields; however, the increase in temperature year by year has a significant negative effect on bean yields; (2) carbon emissions do not directly affect bean production at present but may have an indirect impact on bean production; (3) artificial irrigation and fertilization behavior on bean production has basically reached saturation, making it difficult to continue to increase bean yields and (4) the development of technology and human activity is a mixed blessing, and the consequent inhibiting effects on bean production are currently unable to offset their promoting effects. Thus, when it comes to bean cultivation, China should focus mainly on the overall impact of environmental changes on its production, rather than technical enhancements such as irrigation and fertilization.

**Keywords:** yield per unit area of beans; climate change; panel spatial error model

#### **1. Introduction**

In recent decades, the global climate has undergone significant changes due to natural environmental changes and human activities, and this change is mainly characterized by global warming. The main reason for the warming is the emission of greenhouse gases such as carbon dioxide caused by human activities [1]. Climate change research intends to assess the effects of climate change on ecosystems, economies and societies by predicting possible future changes of climate [2]. The influence of climate change on human beings is comprehensive, multi-level and multi-scale. Agriculture is the sector most directly affected by and vulnerable to climate change. It is also the sector on which human society depends for the output of basic means of living, and the sustainable development of agriculture is directly related to the survival and development of human society [3,4]. As grain is a crucial agricultural product, grain security, grain varieties and grain production have always been the top priority of the global economy. With the advancement of agricultural technologies and the improvement of life quality, the yields of crops and productivity have been boosted over time. At the same time, climate warming, precipitation instability, extreme weather and disasters have brought negative impacts on food production [5]. Studies have shown that the determinants of crop production are climatic (e.g., rainfall, high temperatures) and weather extremes (e.g., floods, droughts and storms) [6]. Changes in rainfall patterns and increased temperature could have a marked impact on food production [7]. According to

**Citation:** Li, S.; You, S.; Song, Z.; Zhang, L.; Liu, Y. Impacts of Climate and Environmental Change on Bean Cultivation in China. *Atmosphere* **2021**, *12*, 1591. https://doi.org/10.3390/ atmos12121591

Academic Editors: Zengyun Hu, Xuguang Tang and Qinchuan Xin

Received: 14 October 2021 Accepted: 26 November 2021 Published: 29 November 2021

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

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

IPCC (2012), crop production is sensitive to climate change associated with temperature increases, changes in rainfall patterns and extreme weather events. In addition to the changes of climate and temperature, carbon dioxide also affects the food production, and it has a higher impact on C3 species (including wheat, rice and bean) than C4 species (including maize and sorghum) [8]. Food production in countries around the world is currently under serious threat from climate change, which leads to challenges to sustainable development as well [9]. Climate change is also an important challenge for food production and security in China. In the context of global warming, changes in temperature and precipitation patterns could lead to a 20–36% decline in maize, wheat, and rice yields in China over the next 20–80 years [10]. China's grain production with the background of climate change has the following three main problems: Firstly, grain production becomes more volatile; secondly, the structure and layout of food production may change, and crop cultivation systems need to change accordingly and thirdly, climate change will substantially increase the cost of grain production [11].