# Long-Term (2007 to 2018) Energy and CO2 Fluxes over an Agriculture Ecosystem in the Southeastern Margin of the Tibetan Plateau

Authors: Anlun Xu; Jian Li; Qun Du; Baoju Dong
Journal: Advances in Meteorology
 (2022)

---

## Abstract

Long-term eddy covariance flux observations over complex topography are crucial for improving the understanding of the turbulent exchanges between the land and atmosphere. Based on a 12-year (2007–2018) record dataset measured with the eddy covariance technique over the Dali agriculture ecosystem in the southeastern margin of the Tibetan Plateau, we investigated the diurnal, seasonal, and interannual variations of the sensible heat flux (Hs), latent heat flux (LE), and carbon dioxide flux (Fc), and their controlling variables. The results showed that Hs and LE exhibited similar diurnal and seasonal variations, while the amplitude of LE was clearly larger than that of Hs throughout the year. The turbulent fluxes showed remarkable fluctuation on the annual scale. The annual average Hs (LE) increases (decreases) from approximately 8 (110) W·m−2 during 2007–2013 to 20 (79) W·m−2 during 2014–2018. The annual cumulative net CO2 ecosystem exchange (NEE) increases significantly from approximately −739 g·C·m−2·yr−1 during 2007–2013 to −218 g·C·m−2·yr−1 during 2014–2018. The relationship between turbulent fluxes and meteorological variables was also examined. Wind speed (WS) is found to be the dominant controlling factor for the Hs on different temporal scales and their correlation coefficients increase when the timescales vary from daily to annual scale; whereas the product of WS and vapor pressure deficit (VPD) is the major meteorological variable controlling the LE over all temporal scales. The net radiation (Rn) is the dominating factor for Fc on daily and monthly timescales, while WS becomes the most controlling factor for Fc on an annual scale. Notably, surface energy and CO2 fluxes are also greatly influenced by the vegetation cover surrounding the measurement site.


---

## Body

## 1. Introduction

The Tibetan Plateau (TP), known as the “Roof of the World,” stretches approximately 2500 km in the longitudinal direction and 1000 km in the latitudinal direction, and its average elevation exceeds 4 km above sea level (ASL), penetrating into the midlevel troposphere (Figure1(a)). Owing to its special atmosphere-land interactions, along with mechanical and thermodynamic effects, the TP exerts significant impacts on the evolution of atmospheric circulation, climate change, and extreme weather events [1–5].Figure 1
(a) Topographic map of the Tibetan Plateau (TP) and Yunnan-Guizhou Plateau (YGP), (b) location of the study area (map from Google Earth), pictures taken on 16 January 2009 (c), 26 September 2010 (d), 10 April 2016 (e), and 17 August 2016 (f) of the Dali flux station.A series of field campaigns and scientific experiments have been carried out over the TP area since the 1970s [6–8] to fully understand the atmospheric scientific questions on the TP. For instance, during May–August in 1979, Chinese scientists performed the Qinghai-Xizang Plateau Meteorological Experiment (QXPMEX) for investigating the radiation and heat budget of the TP as well as its thermal dynamic effects on Asian general circulation [9]. During May–August in 1998, Chinese scientists implemented the Tibetan Plateau Field Experiment (TIPEX) for revealing the plateau’s atmosphere-land physical processes and the structure of the atmospheric boundary layer over the TP [10]. During 1996–2000 and 2001–2005, Chinese, Japanese, and Korean scientists successively carried out the Global Energy and Water Cycle Exchange (GEWEX) Asian Monsoon Experiment (GAME/Tibet), and the Coordinated Enhanced Observing Period (CEOP) Asia-Australia Monsoon Project (CAMP/Tibet) for investigating surface energy exchanges and so on [11, 12]. From 2005 to 2009, Chinese and Japanese scientists carried out the Japan International Cooperation Agency Project (JICA/Tibet Project) for understanding the processes of the land-atmosphere interaction and their effects on the severe weather and climate over the TP and surrounding region [13, 14]. In 2013, the China Meteorological Administration, the National Natural Science Foundation of China, and the Chinese Academy of Sciences jointly initiated the third Tibetan Plateau Atmospheric Scientific Experiment (TIPEX-III) and formally launched in 2014 for understanding the land-surface boundary layer and cloud-precipitation physical processes and tropospheric, stratospheric, and atmospheric composition exchanges over the TP, with an 8–10-year implementation plan [8, 15]. The surface energy budget, turbulent exchanges, and atmospheric boundary layer were regarded as the important observational objectives and research tasks during these field campaigns. Based on the in situ observation datasets, many studies have reported the characteristics of surface water vapor, heat energy, and CO2 fluxes.The surface over the TP is an atmospheric heat source throughout the year [16], with the heat source being mainly regulated by the sensible heat in the dry season (November to April of the next year), where as it being modulated by the release of latent heat in the wet season (May to October) [17–19]. The sensible heat flux (Hs) over the TP is dominant under the south branch of the westerly wind field, while the latent heat flux (LE) is dominant under the plateau monsoon circulation field [20]. The Hs gradually decreases from western TP to the Chengdu plain, and LE increases significantly from the western to central TP [21]. The daytime dominance of Hs and LE in the eastern edge of the TP varies seasonally with LE dominating in summer (June–August) and Hs dominating in winter (December–February) [22]. The energy imbalance over the TP was obvious with the surface energy balance ratio of 0.74 in summer [23]. In addition, lots of long-term (over 10 years) fluxes observational stations, such as those in alpine shrubland [24], steppe [25], meadow, and desert [26], have already been established on the TP with the support of the abovementioned experiments or other projects.The southeastern TP belongs to the transitional zone between the TP and Yunnan-Guizhou Plateau (YGP), which is located in a major water vapor path of southwest China and a confluence zone between the South Asian monsoon and East Asian monsoon [27], and where it is also sensitive to atmospheric adiabatic heating over the TP. This area has a complex topography with mountains, lakes, rivers, basins, meadows, forests, farmlands, and wetlands. Due to the unique location, complex terrain, and a lack of in situ observation data in this area, the output products of some atmospheric general circulation models show obvious deviations [28]. Moreover, basic information on the land surface and atmosphere in this region is crucial for the evolution of Asian atmospheric circulation as well as weather and climate change in the Yangtze River basin and southwestern China [29, 30]. Therefore, revealing the characteristics of the local weather and climate in this region is of great importance, and the quantification of atmosphere-land interaction parameters is helpful for improving the numerical models and parameterization schemes of the atmospheric boundary layer in southwest China and even in the East Asia region.Knowledge on the energy and carbon exchange process in the southeastern TP is still very limited due to its complicated topography. In order to improve the atmosphere-land interaction in this region, a long-term observation of energy and carbon fluxes were conducted in a cropland in the southeastern border of the TP. The flux measurements were based on the eddy covariance technique. This observation site (100.177°E, 25.709°N, 1977.7 m ASL, Figures 1(b)–1(f)), which was built with the support of the JICA/Tibet Project in December 2006 [13, 14], is in a basin between Cangshan Mountain (CM, also known as “Diancangshang Mountain”) and Yu’anshan Mountain. The Erhai Lake also lies in the basin. The Erhai Lake and CM are about 2 km and 4 km away from the cropland observation site, respectively. The long-term observation provides a unique opportunity to investigate the patterns of energy and carbon fluxes, which will be of great help to evaluate and improve the atmosphere-land interaction parameters over the southeastern margin of the TP. The objectives of this study are to examine the diurnal, seasonal, and interannual variability in surface energy and CO2 fluxes and determine their meteorological controlling factors in different temporal scales over the southeast extension of the TP.


## 2. Materials and Methods

### 2.1. Station

The planetary boundary layer (PBL) flux station, which is located in Dali Bai Autonomous Prefecture, Yunnan Province, China, is surrounded by an open and flat farmland area between the CM and Erhai Lake (Figure1(b)) where the soil is principally sandy loam. During the observational period, this cropland was planted with broad beans and rice in rotation until November 2013, that is, broad beans were mainly planted in the dry season and rice was mainly planted in the wet season from January 2007 to October 2013 (Figures 1(c)–1(d)). The maximum broad bean height at peak growth can reach up to 1.0 m, and the rice height can reach up to 1.2 m. After November 2013, small trees were planted on the east and north side of the PBL flux tower croplands (Figures 1(e)–1(f)), with the height of the trees approximately ranging from 1 to 5 m.For assessing the source area of the measured fluxes, the footprint analysis is conducted by using the Kormann and Meixner model [31].(1)fyx,z=1Γμξμx1+μe−ξ/x,μ=1+mr,ξ=ωzrr2q,where Γ, μ, and ξ are the gamma function, a constant, and flux length scale, respectively. ω and q are a constant in the wind speed power law profile and in the turbulent diffusion power law profile, respectively. r is a shape factor (where r=m+2−n, and m and n are indexes of the wind speed and turbulent diffusion power law profiles, respectively). According to the Kormann and Meixner model analysis, it can be found that 95% of the source area contributing to the measured fluxes during 2007–2018 was approximately shaped like an ellipse, which mainly came from 700 and 500 m in the southeast and northwest directions of the PBL flux tower (Figure 2), relying on the predominant wind direction. Although several roads, trees, buildings, and irrigation channels were located in the source area, the underlying surface of the source area was primarily covered by croplands, which satisfied the requirement of eddy covariance measurements. It indicated that the measured fluxes were reliable and mainly contributed by the croplands, and the measured fluxes had good representative.Figure 2
Average footprint source at Dali flux station from 2007 to 2018. The cumulative footprint contours were superimposed on the map from Google Earth and the maximum footprint contour line is shown with 95%.Dali generally experiences a windy and dry season influenced by the southward movement of the westerlies, and a wet season controlled by the warm and moist airflow coming from the Bay of Bengal in the southwest direction and from the subtropical Pacific high in the southeast direction [32]. On average, the wet season in the study area lasts from late May to the end of early October every year.


### 2.2. Measurement Setup

The instrumentations, measurement heights, sampling rate, sensors, and manufacturers of the Dali flux station is shown in Table1. A 20 m tower was set up to acquire the vertical profiles of air temperature (Ta), relative humidity (RH), wind speed (WS), and wind direction (WD) along with the turbulent fluxes of Hs, LE, and CO2 fluxes (Fc) within the near-surface layer. The Ta and RH at four heights were measured with four temperature and relative humidity probes, and the WS and WD at the same heights were measured with four wind sensors. An eddy covariance (EC) system was placed at 5.08 m height on the tower to directly determine the Hs, LE, and Fc, which contained a three-dimensional sonic anemometer and an open-path CO2/H2O infrared gas analyser. A net radiometer was used to measure the net radiation flux (Rn) as well as downward/upward longwave/shortwave radiation components. The downward shortwave radiation was taken as the global solar radiation (Rg) in this study. An infrared temperature sensor was used to measure the soil surface temperature. A tipping bucket rain gauge was used to measure the precipitation (PPT). Five temperature probes and water content reflectometers were used to measure the soil temperature (Ts) and soil water content (SWC) at the same depths, respectively. Three soil heat flux plates were used to measure the soil heat flux (Gs). A digital barometer was used to measure the barometric pressure. The data of three-dimensional wind speed, sonic temperature, CO2 density, and H2O density were recorded by a data logger with a 1 GB CF card. The rest data were recorded by the same data logger on a 64 MB CF card. More observational details about this flux station were documented in our previous study [33].Table 1
Instrumentations, measurement heights, sampling rate, sensors, and manufacturers of the Dali flux station.
InstrumentsMeasurement heights (m)Sampling rateSensors and manufacturersThree-dimensional sonic anemometer5.0810 HzCSAT3, Campbell Scientific Inc., Logan, UT, USAOpen-path CO2/H2O infrared gas analyser5.0810 HzLI-7500, LI-COR Inc., Lincoln, NE, USATemperature and relative humidity probe2, 4, 10, 2010 minHMP45C, Vaisala, Helsinki, FinlandWind sensor2, 4, 10, 2010 min034B, Met One Instruments Inc., Grants Pass, OH, USANet radiometer1.510 minCNR1, Kipp & Zonen, Delft, the NetherlandsInfrared temperature sensor1.510 minIRR-P, Apogee Instruments Inc., Logan, UT, USABucket rain gauge0.310 minRG13H, Vaisala, Helsinki, FinlandSoil temperature probe−0.04, −0.1, −0.2, −0.6, −1.010 minModel 107, Campbell Scientific Inc., Logan, UT, USAWater content reflectometer−0.04, −0.1, −0.2, −0.6, −1.010 minCS616, Campbell Scientific Inc., Logan, UT, USASoil heat flux plate−0.04, −0.1, −0.210 minHFP01, Hukseflux Thermal Sensors, Delft, the NetherlandsDigital barometer1.010 minPTB220, Vaisala, Helsinki, FinlandData logger1.0CR3000, Campbell Scientific Inc., Logan, UT, USA


### 2.3. Data Processing

To obtain 30 min eddy covariance flux data (Hs, LE, and Fc), the EddyPro software (version 6.2.1, 2017, LI-COR, Inc., Lincoln, NE, USA) was applied to process the 10 Hz data. All processing settings were implemented in the EddyPro software (LI-COR, Inc.). The raw data processing and flux correcting steps included spike detection/removal [34], axis rotation (double rotation method) for tilt correction [35], high and low frequency spectral correction [36, 37], and Webb–Pearman–Leuning (WPL) correction [38]. The turbulent fluxes were calculated as follows:(2)Hs=ρCpw′T′¯,LE=ρLvw′q′¯,Fc=ρw′c′¯,where ρ, Cp, and Lv are the air density (kg·m−3), the heat capacity of air at constant pressure (J·kg−1·K−1), and the latent heat of vaporization of water (J·kg−1), respectively. w′, T′, q′, and c′ are the fluctuations of the vertical wind component (m·s−1), air temperature (K), specific humidity (g·kg−1), and the concentration of CO2 (μmol·m−3), respectively.In addition, a quality check of the 30 min fluxes was performed using the steps proposed by Foken et al., including the steady state and well-developed turbulent conditions test [39]. The calculated flux values with low quality flag as well as instrument malfunction and under rainy conditions were discarded. After data quality check, the high and moderate quality data for Hs, LE, and Fc accounted for 76.4%, 70.4%, and 71.3%, respectively. Based on the REddyPro R package (Max Planck Institute for Biogeochemistry, https://www.bgc-jena.mpg.de/bgi/index.php/Services/REddyProcWebRPackage), the gap filling of flux data and meteorological data was performed in the cross-platform language R [40], using methods similar to Falge et al. [41] and adding the covariance between flux and meteorological data as well as the temporal autocorrelation of the flux [42]. The percentage of Hs, LE, and Fc data that required gap-filling are given in Table 2. Ancillary meteorological data included the friction velocity (u∗), Rg, Ts at 4 cm soil depth, Ta, and RH at 4 m height. The vapor pressure deficit (VPD) was calculated with Ta and RH. The times given in our study were in Beijing Time (UTC+8).Table 2
Percentage of turbulent flux data that required gap-filling.
Hs (%)LE (%)Fc (%)11.618.318.6To access the EC data quality, the surface energy balance ratio (EBR) was driven from the data of turbulent energy flux (Hs + LE) and available energy flux (Rn − Gs), using the linear regression method.Based on the daily and monthly averaged data from 2007 to 2018, the relationships between Hs + LE and Rn − Gs are shown in Figure3. The daily and monthly averaged EBR are 0.81 and 0.79 over Dali cropland during the whole measurement period, respectively. It is reported that the EBR in majority flux sites ranged from 70% to 90% [43]. The EBR value in this study was comparable to the FLUXNET sites. The surface energy imbalance was caused by many reasons, for example, the instrument errors (such as systematic error of radiation measurements and mismatch of footprints), data processing errors (such as uncertainty due to averaging operator), additional sources of energy (such as canopy and soil storage), submesoscale transport processes (such as vertical/horizontal advection) [44, 45].Figure 3
Relationship between the daily (a) and monthly (b) averaged turbulent energy flux (Hs+ LE) and available energy flux (Rn − Gs) from 2007 to 2018.
(a)(b)The normalized difference vegetation index (NDVI) data during 2007–2018 was obtained from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1) products of the National Aeronautics and Space Administration (NASA) [46]. The temporal and spatial resolution of these products are 16-day and 250 m, respectively. The averaged NDVI value of the four points surrounding the measurement location was regarded as the vegetation cover in this flux station. Due to the lack of biometric measurements (e.g., leaf area index) and other detailed records, NDVI was applied to analyse the development of vegetation.To quantify the controlling factor of turbulent fluxes, correlation analysis of the Hs, LE, and Fc with environmental variables was undertaken using the IBM SPSS statistics software (version 25.0, IBM Corp., Chicago, IL, USA). Furthermore, the topographic map and figures were drawn using ArcGIS software (version 10.4.1, Esri Inc., Redlands, CA, USA) and OriginPro software (version 2018, OriginLab Corp., Northampton, MA, USA), respectively.


## 2.1. Station

The planetary boundary layer (PBL) flux station, which is located in Dali Bai Autonomous Prefecture, Yunnan Province, China, is surrounded by an open and flat farmland area between the CM and Erhai Lake (Figure1(b)) where the soil is principally sandy loam. During the observational period, this cropland was planted with broad beans and rice in rotation until November 2013, that is, broad beans were mainly planted in the dry season and rice was mainly planted in the wet season from January 2007 to October 2013 (Figures 1(c)–1(d)). The maximum broad bean height at peak growth can reach up to 1.0 m, and the rice height can reach up to 1.2 m. After November 2013, small trees were planted on the east and north side of the PBL flux tower croplands (Figures 1(e)–1(f)), with the height of the trees approximately ranging from 1 to 5 m.For assessing the source area of the measured fluxes, the footprint analysis is conducted by using the Kormann and Meixner model [31].(1)fyx,z=1Γμξμx1+μe−ξ/x,μ=1+mr,ξ=ωzrr2q,where Γ, μ, and ξ are the gamma function, a constant, and flux length scale, respectively. ω and q are a constant in the wind speed power law profile and in the turbulent diffusion power law profile, respectively. r is a shape factor (where r=m+2−n, and m and n are indexes of the wind speed and turbulent diffusion power law profiles, respectively). According to the Kormann and Meixner model analysis, it can be found that 95% of the source area contributing to the measured fluxes during 2007–2018 was approximately shaped like an ellipse, which mainly came from 700 and 500 m in the southeast and northwest directions of the PBL flux tower (Figure 2), relying on the predominant wind direction. Although several roads, trees, buildings, and irrigation channels were located in the source area, the underlying surface of the source area was primarily covered by croplands, which satisfied the requirement of eddy covariance measurements. It indicated that the measured fluxes were reliable and mainly contributed by the croplands, and the measured fluxes had good representative.Figure 2
Average footprint source at Dali flux station from 2007 to 2018. The cumulative footprint contours were superimposed on the map from Google Earth and the maximum footprint contour line is shown with 95%.Dali generally experiences a windy and dry season influenced by the southward movement of the westerlies, and a wet season controlled by the warm and moist airflow coming from the Bay of Bengal in the southwest direction and from the subtropical Pacific high in the southeast direction [32]. On average, the wet season in the study area lasts from late May to the end of early October every year.


## 2.2. Measurement Setup

The instrumentations, measurement heights, sampling rate, sensors, and manufacturers of the Dali flux station is shown in Table1. A 20 m tower was set up to acquire the vertical profiles of air temperature (Ta), relative humidity (RH), wind speed (WS), and wind direction (WD) along with the turbulent fluxes of Hs, LE, and CO2 fluxes (Fc) within the near-surface layer. The Ta and RH at four heights were measured with four temperature and relative humidity probes, and the WS and WD at the same heights were measured with four wind sensors. An eddy covariance (EC) system was placed at 5.08 m height on the tower to directly determine the Hs, LE, and Fc, which contained a three-dimensional sonic anemometer and an open-path CO2/H2O infrared gas analyser. A net radiometer was used to measure the net radiation flux (Rn) as well as downward/upward longwave/shortwave radiation components. The downward shortwave radiation was taken as the global solar radiation (Rg) in this study. An infrared temperature sensor was used to measure the soil surface temperature. A tipping bucket rain gauge was used to measure the precipitation (PPT). Five temperature probes and water content reflectometers were used to measure the soil temperature (Ts) and soil water content (SWC) at the same depths, respectively. Three soil heat flux plates were used to measure the soil heat flux (Gs). A digital barometer was used to measure the barometric pressure. The data of three-dimensional wind speed, sonic temperature, CO2 density, and H2O density were recorded by a data logger with a 1 GB CF card. The rest data were recorded by the same data logger on a 64 MB CF card. More observational details about this flux station were documented in our previous study [33].Table 1
Instrumentations, measurement heights, sampling rate, sensors, and manufacturers of the Dali flux station.
InstrumentsMeasurement heights (m)Sampling rateSensors and manufacturersThree-dimensional sonic anemometer5.0810 HzCSAT3, Campbell Scientific Inc., Logan, UT, USAOpen-path CO2/H2O infrared gas analyser5.0810 HzLI-7500, LI-COR Inc., Lincoln, NE, USATemperature and relative humidity probe2, 4, 10, 2010 minHMP45C, Vaisala, Helsinki, FinlandWind sensor2, 4, 10, 2010 min034B, Met One Instruments Inc., Grants Pass, OH, USANet radiometer1.510 minCNR1, Kipp & Zonen, Delft, the NetherlandsInfrared temperature sensor1.510 minIRR-P, Apogee Instruments Inc., Logan, UT, USABucket rain gauge0.310 minRG13H, Vaisala, Helsinki, FinlandSoil temperature probe−0.04, −0.1, −0.2, −0.6, −1.010 minModel 107, Campbell Scientific Inc., Logan, UT, USAWater content reflectometer−0.04, −0.1, −0.2, −0.6, −1.010 minCS616, Campbell Scientific Inc., Logan, UT, USASoil heat flux plate−0.04, −0.1, −0.210 minHFP01, Hukseflux Thermal Sensors, Delft, the NetherlandsDigital barometer1.010 minPTB220, Vaisala, Helsinki, FinlandData logger1.0CR3000, Campbell Scientific Inc., Logan, UT, USA


## 2.3. Data Processing

To obtain 30 min eddy covariance flux data (Hs, LE, and Fc), the EddyPro software (version 6.2.1, 2017, LI-COR, Inc., Lincoln, NE, USA) was applied to process the 10 Hz data. All processing settings were implemented in the EddyPro software (LI-COR, Inc.). The raw data processing and flux correcting steps included spike detection/removal [34], axis rotation (double rotation method) for tilt correction [35], high and low frequency spectral correction [36, 37], and Webb–Pearman–Leuning (WPL) correction [38]. The turbulent fluxes were calculated as follows:(2)Hs=ρCpw′T′¯,LE=ρLvw′q′¯,Fc=ρw′c′¯,where ρ, Cp, and Lv are the air density (kg·m−3), the heat capacity of air at constant pressure (J·kg−1·K−1), and the latent heat of vaporization of water (J·kg−1), respectively. w′, T′, q′, and c′ are the fluctuations of the vertical wind component (m·s−1), air temperature (K), specific humidity (g·kg−1), and the concentration of CO2 (μmol·m−3), respectively.In addition, a quality check of the 30 min fluxes was performed using the steps proposed by Foken et al., including the steady state and well-developed turbulent conditions test [39]. The calculated flux values with low quality flag as well as instrument malfunction and under rainy conditions were discarded. After data quality check, the high and moderate quality data for Hs, LE, and Fc accounted for 76.4%, 70.4%, and 71.3%, respectively. Based on the REddyPro R package (Max Planck Institute for Biogeochemistry, https://www.bgc-jena.mpg.de/bgi/index.php/Services/REddyProcWebRPackage), the gap filling of flux data and meteorological data was performed in the cross-platform language R [40], using methods similar to Falge et al. [41] and adding the covariance between flux and meteorological data as well as the temporal autocorrelation of the flux [42]. The percentage of Hs, LE, and Fc data that required gap-filling are given in Table 2. Ancillary meteorological data included the friction velocity (u∗), Rg, Ts at 4 cm soil depth, Ta, and RH at 4 m height. The vapor pressure deficit (VPD) was calculated with Ta and RH. The times given in our study were in Beijing Time (UTC+8).Table 2
Percentage of turbulent flux data that required gap-filling.
Hs (%)LE (%)Fc (%)11.618.318.6To access the EC data quality, the surface energy balance ratio (EBR) was driven from the data of turbulent energy flux (Hs + LE) and available energy flux (Rn − Gs), using the linear regression method.Based on the daily and monthly averaged data from 2007 to 2018, the relationships between Hs + LE and Rn − Gs are shown in Figure3. The daily and monthly averaged EBR are 0.81 and 0.79 over Dali cropland during the whole measurement period, respectively. It is reported that the EBR in majority flux sites ranged from 70% to 90% [43]. The EBR value in this study was comparable to the FLUXNET sites. The surface energy imbalance was caused by many reasons, for example, the instrument errors (such as systematic error of radiation measurements and mismatch of footprints), data processing errors (such as uncertainty due to averaging operator), additional sources of energy (such as canopy and soil storage), submesoscale transport processes (such as vertical/horizontal advection) [44, 45].Figure 3
Relationship between the daily (a) and monthly (b) averaged turbulent energy flux (Hs+ LE) and available energy flux (Rn − Gs) from 2007 to 2018.
(a)(b)The normalized difference vegetation index (NDVI) data during 2007–2018 was obtained from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1) products of the National Aeronautics and Space Administration (NASA) [46]. The temporal and spatial resolution of these products are 16-day and 250 m, respectively. The averaged NDVI value of the four points surrounding the measurement location was regarded as the vegetation cover in this flux station. Due to the lack of biometric measurements (e.g., leaf area index) and other detailed records, NDVI was applied to analyse the development of vegetation.To quantify the controlling factor of turbulent fluxes, correlation analysis of the Hs, LE, and Fc with environmental variables was undertaken using the IBM SPSS statistics software (version 25.0, IBM Corp., Chicago, IL, USA). Furthermore, the topographic map and figures were drawn using ArcGIS software (version 10.4.1, Esri Inc., Redlands, CA, USA) and OriginPro software (version 2018, OriginLab Corp., Northampton, MA, USA), respectively.


## 3. Results and Discussion

### 3.1. Meteorological Conditions and Vegetation Cover

The climate at the measurement site features a monsoon climate of the low latitude plateau. Based on the observational data at an automatic weather station (AWS) located approximately 220 m southwest of the study site, the monthly average Ta during 1951–2018 ranged from 8.6°C in January to 20.2°C in June, with an annual average Ta of 15.1°C (Figure4). The monthly average PPT during 1951–2018 ranged from 12.9 mm in December to 218.4 mm in August, with the annual average PPT of 1048.1 mm and nearly 85% of precipitation falling in May to October. The monthly average WS during 1951–2018 ranged from 1.4 m·s−1 in September to 3.7 m·s−1 in March with the annual average WS of 2.4 m·s−1 and the strongest WS of 40.8 m·s−1. The gale is frequently experienced from winter through spring.Figure 4
Box chart of monthly averaged air temperature (a) and precipitation (b) from 1951 to 2018. The circle, square, and five-pointed star represent the mean, minimum, and maximum values, respectively. The horizontal line on the box represents the median and the box represents the 25–75th percentiles. The whisker represents the 5–95th percentile.
(a)(b)The meteorological variables during the study period are given in Table3 and shown in Figure 5. The difference of every variable in the dry and wet season is significant, namely, Rg, Ta, and RH, in the wet season are higher than those in the dry season. The solar energy resources are rich, with the annual integrated Rg varying from 5066 to 6703 MJ·m−2·yr−1, and a 12-year average value of 6154.1 MJ·m−2. The daily average Ta at heights of 2, 4, 10, and 20 m are always above 0°C, with the maximum and minimum values recorded in June 2014 and December 2013. The daily average Ts at the upper soil layer are higher than those at the lower soil layer in the wet season; conversely, Ts at the lower soil layer are higher than those at the upper soil layer in the dry season. The monthly average Ts values reach the maximum in July and the lowest in January. The RH within the near-surface layer usually decreases with increasing measurement height. The monthly average RH values reach the highest in August and the lowest in February. The WS within the near-surface layer commonly increases with increasing measurement height. The monthly average WS values reach the strongest in February and the weakest in September. The SWC values respond significantly to PPT. With the beginning of the wet season, the SWC rapidly increases and reaches its maximum.Table 3
Annual total global solar radiation (Rg), annual averaged wind speed (WS) at 10 m height, air temperature (Ta), and relative humidity (RH) at 4 m height, vapor pressure deficit (VPD), soil temperature (Ts), and soil water content (SWC) at 4 cm depth, and annual total precipitation (PPT) at the Dali site from 2007 to 2018.
YearRg (MJ·m−2)WS (m·s−1)Ta (°C)RH (%)VPD (kPa)Ts (°C)SWC (m−3·m−3)PPT (mm)20075066.12.616.668.00.6618.00.3201139.820085659.32.515.868.30.6218.40.3811322.820096702.72.716.462.10.7617.30.3521031.220106043.62.816.362.60.7416.40.394863.820116126.72.515.762.40.7215.90.426740.120126120.72.716.457.70.8315.80.422853.220136154.72.415.957.60.8416.10.423732.920146687.12.416.558.00.8515.90.332812.920156440.72.016.161.00.8115.90.326984.920166069.81.915.964.80.7015.30.3741153.620176356.92.016.163.60.7316.40.381909.820186421.61.915.563.00.7216.90.3441000.7Average6154.12.416.162.40.7516.50.373962.1Standard deviation450.20.30.33.60.070.90.038179.3Figure 5
Daily integrated global solar radiation (Rg) (a), daily averaged air temperature (Ta) (b), relative humidity (RH) (c), and wind speed (WS) (d) at four heights, daily averaged soil temperature (Ts) at five depths (e), and daily averaged soil water content (SWC) at five depths and daily total precipitation (PPT) (f) from 2007 to 2018.
(a)(b)(c)(d)(e)(f)The maximum daily total PPT value is observed in October 2015 (121.0 mm·d−1), and the maximum monthly total PPT is observed in June 2008 (336.7 mm·month−1). The annual total PPT fluctuates from 732.9 mm·yr−1 (2013) to 1322.8 mm·yr−1 (2008), with an average value of 962.1 mm (Table 3). The PPT in the wet season ranges from 643.2 to 1171.9 mm, accounting for 75.9% to 93.6% of the annual total PPT. Because of the frequent drought influence, the PPT for five consecutive years from 2010 to 2014 is below the 12-year average value.Local complicated terrain often affects the airflow within the atmospheric boundary layer, which not only changes wind speed but also changes wind direction. The predominant wind direction is easterly and east-southeasterly throughout the year (FigureS1). The prevailing and subprevailing winds at 10 m height display a clear diurnal cycle, with east-southeasterly and easterly winds during the daytime and west-northwesterly and static winds at night (Figure S2). There is a good correlation between the time when the prevailing winds switches and the sunrise and sunset (solar radiation heating effect). This result may be chiefly caused by the orientation of the CM and Erhai Lake, which extends from northwest to southeast.The NDVI has an obvious change before and after November 2013 because of the surface alternation (Figure6). For each year, the NDVI has two peaks in January or February (the fastest growing period of broad bean) and in July or August (the fastest growing period of rice), and it has also two valleys in April or May and in September or October (the harvesting periods of broad bean and rice) from 2007 to 2013. However, the NDVI has only a peak in July or August (growing period of the mixing plants) and a valley in January or February (nongrowing period of the mixing plants) from 2014 to 2018.Figure 6
Seasonal and interannual variations of the normalized difference vegetation index (NDVI) at 250 m spatial resolution and 16-day temporal resolution from 2007 to 2018.


### 3.2. Surface Energy Fluxes

#### 3.2.1. Diurnal, Seasonal, and Interannual Variations in Surface Energy Fluxes

The monthly average diurnal variations in the energy balance components from 2007 to 2018 are shown in Figures7(a)–7(d). As expected, both heat and evapotranspiration vary with available net radiation. For each month, the Rn, Hs, LE, and Gs all exhibit clear, similar diurnal courses. During the day, Rn mainly comes from solar shortwave radiation and remains positive values, reaching its peak when the sun is in the zenith (between 12 : 00 and 13 : 00, depending on the season), while Rn chiefly comes from surface longwave radiation and becomes a negative value at night, reaching its minimum between 19 : 00 and 21 : 00. The monthly valleys and peaks of hourly average Rn vary from −67 W·m−2 in February to −32 W·m−2 in October and from 418 W·m−2 in December to 513 W·m−2 in June, respectively. Hs also attains its maximum and minimum when Rn reaches the peak and valley, which is consistent with the results found by Acharya et al. [47] for an agricultural field in Nepal. The monthly maximal and minimal Hs values vary from 66 W·m−2 in July to 127 W·m−2 in February and from −52 W·m−2 in March to −7 W·m−2 in August. The peak times of both LE and Gs are approximately 1 hour later than those of Rn and Hs. Similarly, Ajao et al. [48] reported that there was a lag of about 1 hour between Hs and LE at an agricultural site in Nigeria due to a phase change (water to vapor). LE remains a positive value and implies evaporation throughout the day. The monthly maximal and minimal LE values vary from 205 W·m−2 in December to 337 W·m−2 in August and from 3 W·m−2 in January to 20 W·m−2 in July. The diurnal average Gs is relatively small throughout the year, with its valley occurring in November (−25 W·m−2) and the peak occurring in June (47 W·m−2). The diurnal variations of Hs and LE in this study is obviously different from that observed in Erhai Lake [49], which may be related to the difference between land and water surface.Figure 7
Monthly average diurnal variations in net radiation flux (Rn) (a), sensible heat flux (Hs) (b), latent heat flux (LE) (c), and soil heat flux (Gs) (d), and monthly average diurnal variation in CO2 flux (Fc) (e) from 2007 to 2018. The error bars indicate plus or minus one standard deviation.
(a)(b)(c)(d)(e)The seasonal and interannual variations in the monthly average energy balance components from 2007 to 2018 are shown in Figures8(a)–8(d) and 9, respectively. For each year, the Rn, LE, and Gs values in the wet season are all larger than those in the dry season. The monthly average Rn generally attains the peak in June or July and the valley in January or December, with the maximum and minimum values of 229.0 and 22.8 W·m−2. The monthly average Hs shows a significant seasonal variation from 2007 to 2013, while Hs remains high and does not present evident variability from 2014 to 2018. Meanwhile, the annual average Hs increases from approximately 6 W·m−2 during 2007–2012 to 19 W·m−2 during 2013–2018. Although LE shows strong changes between dry and wet seasons during the entire period, the seasonal change amplitude of LE becomes small after 2014, and the annual average LE decreases from approximately 110 W·m−2 during 2007–2013 to 79 W·m−2 during 2014–2018. Overall, there is an obvious step change for Hs after 2013 and for LE after 2014, which may mainly be regulated by the vegetation types and cover around the PBL flux tower.Figure 8
Seasonal and interannual variations in monthly averaged net radiation flux (Rn) (a), sensible heat flux (Hs) (b), latent heat flux (LE) (c), soil heat flux (Gs) (d), and CO2 flux (Fc) (e) from 2007 to 2018.
(a)(b)(c)(d)(e)Figure 9
Annual averaged sensible heat flux (Hs), latent heat flux (LE), net radiation flux (Rn), and soil heat flux (Gs) from 2007 to 2018.


#### 3.2.2. Surface Energy Partition

Owing to the change of vegetation cover, the proportions of the sensible and latent heat fluxes in net radiation (Hs/Rn and LE/Rn) are distinctly different. The annual average Hs/Rn increased from 5.3% during 2007–2012 to 17.7% during 2013–2018, while the annual average LE/Rn decreased from 94.6% during 2007–2013 to 72.8% during 2014–2018. In addition, the annual average Hs/Rn and LE/Rn during 2007–2018 are 11.2% and 85.8%, respectively. The annual average Bowen ratio (Hs/LE) during the study period is 0.15, suggesting that the surface heat exchange is modulated by the LE in the whole year. These findings are consistent with the study of Wang et al. [50] who reported that LE was the main part of energy consumption, with its annual value being more than twice of Hs in the karst region of the YGP, but not consistent with Li et al. [7] who found that the major energy source providing heat to the atmosphere was the Hs during the premonsoon period but the LE during the monsoon period at four flux stations (Linzhi, Nagqu, NamCo, and Qomolangma) over the TP. Through this comparative analysis, it is implied that the characteristics of energy exchange over the TP, YGP, and their transitional zones are similar in the wet season but clearly different in the dry season due to the influence of monsoon [7], the difference of location, and land-cover type.Furthermore, the monthly average Gs usually reaches the peak in May or June and the valley in November or December (Figure7(d)). The annual average proportion of the soil heat flux in net radiation (Gs/Rn) decreased from −0.09% during 2007–2013 to −6.3% during 2014–2018.


#### 3.2.3. Controlling Variables in Hs at Different Temporal Scales

The energy, H2O and CO2 exchanges between the atmosphere and croplands depend not only on biological control but also on environmental factors such as meteorological conditions, vegetation characteristics, soil properties, and management practices [51–54]. Relationships between the Hs and environmental variables on timescales from half-hourly to yearly are given in Table 4. Compared to other environmental variables, the WS is the most important controlling variable for Hs from daily to yearly timescales. The correlation coefficients between them become higher as the timescales become longer, with r ranging from −0.95 to −0.55. The temperature difference between the upper soil layer (at 4 cm depth) and the atmosphere (ΔT) as well as WS multiplied by ΔT present positive effects on Hs on daily and monthly timescales but negative effects on Hs at half-hourly and yearly timescales. The correlation coefficients between Hs and WS are higher than those between Hs and ΔT as well as WS multiplied by ΔT, especially on the yearly timescale. Rn has significantly positive effects on Hs from half-hourly to monthly timescales (r ranging from 0.26 to 0.83) but a negative effect at yearly timescale (r = −0.39). The correlation coefficients between Hs and the rest of the environmental factors are relatively low.Table 4
Correlation coefficients of the averaged sensible heat flux (Hs) and environmental variables including the averaged wind speed (WS) at 4 m height, temperature difference between the upper soil layer (at 4 cm depth) and the atmosphere (ΔT), the WS multiplied by ΔT (WS∗ ΔT), net radiation (Rn) at 1.5 m height, air temperature (Ta), the total precipitation (PPT), and the normalized difference vegetation index (NDVI) at multiple timescales from 2007 to 2018. There is no examination about the PPT at half-hourly timescales and the NDVI between half-hourly and daily timescales.
TimescalesEnvironmental variablesWSΔTWS∗ ΔTRnTaPPTNDVIHalf-hourly0.01∗−0.20∗−0.10∗0.83∗0.31∗Daily−0.55∗0.19∗0.28∗0.53∗0.17∗−0.13∗Monthly−0.75∗0.38∗0.48∗0.26∗0.27∗0.25∗−0.08Yearly−0.95∗−0.39−0.32−0.39−0.49−0.100.03Note: the asterisk (∗) indicates the significance level of 0.05.There may be several reasons for the very high correlation of Hs with WS at the annual timescale. On the one hand, the alteration of the vegetation cover results in the surface roughness change, which is high enough to affect significantly the WS, that is, the average WS of 2.0 m·s−1 during 2014 to 2018 is significantly lower than that of 2.6 m·s−1 during 2007 to 2013 (Table 3). At the same time, the Hs significantly increases (Figure 8(b)). The close correlation between WS and Hs was also reported in many other studies. On the other hand, it is well known that the Hs is primarily determined by ΔT as well as by the turbulent exchange coefficient [55]. There is a strong local circulation within the boundary layer existing in the study area, which is caused by the thermally generated airflow (lake-land breeze and mountain-valley breeze), and the dynamically generated airflow (gorge wind) [56, 57]. This local circulation could modify wind speed and direction [33] and result in the air column to be vertically well mixed [58]. It may create more turbulent eddies and promote the vertical exchange of heat, H2O, and CO2 between the atmosphere and cropland when the WS increases but is less than a threshold. Wind-driven mixing largely modifies meteorological factors, which results in the changes of heat and water vapor exchange [59].


#### 3.2.4. Controlling Variables in LE at Different Temporal Scales

The correlation coefficients of the LE and its environmental controls at multiple timescales are also summarized in Table5. Relative to other environmental variables, the product of WS and VPD is the main environmental variable controlling LE, showing a significantly positive correlation over all timescales, with r ranging from 0.41 to 0.72. Although WS and VPD are positively correlated with LE from half-hourly to monthly timescales (r being more than 0.32), the values of r between WS and LE are slightly smaller than those between VPD and LE. On a yearly timescale, WS has a strong positive correlation with LE (r = 0.95), while VPD has a weak negative correlation with LE (r = −0.10). This implies that the effect of WS is not significant as VPD on LE when the timescale is below the monthly timescale, whereas WS plays a decisive role in affecting LE on a yearly scale. Clearly, the impacts of WS and VPD, and their products on LE relies on timescales [60]. There are positive correlations between Rn and LE on all timescales with r ranging from 0.31 to 0.82, which indicates that Rn primarily regulates the variation in water vapor flux, particularly on half-hourly timescales. Ta displays significant positive correlations with LE below monthly timescales but nonsignificant correlation at the yearly timescale.Table 5
Similar to Table4, the flux variable is the latent heat flux (LE), and the environmental variables include the WS, Rn, Ta, PPT, NDVI, as well as the vapor pressure deficit (VPD), and the WS multiplied by VPD (WS ∗ VPD) at 4 m height.
TimescalesEnvironmental variablesWSVPDWS∗VPDRnTaPPTNDVIHalf-hourly0.32∗0.53∗0.43∗0.82∗0.49∗Daily0.32∗0.45∗0.41∗0.66∗0.56∗−0.11∗Monthly0.42∗0.43∗0.45∗0.58∗0.61∗0.21∗0.11Yearly0.95∗−0.100.72∗0.310.34−0.08−0.12


### 3.3. CO2 Flux

#### 3.3.1. Diurnal, Seasonal, and Interannual Variations in CO2 Flux

The monthly average diurnal variations in Fc from 2007 to 2018 are shown in Figure7(e). For each month, the Fc has a noticeable diurnal cycle with positive values at night reaching their maximum at approximately 07 : 00 and negative values during the day reaching their minimum at approximately 13 : 00. This result indicates that the study area acts as a weak carbon source at night because of ecosystem respiration, but acts as a carbon sink during the day as the CO2 uptake sourced from plant photosynthesis significantly exceeds CO2 emission released from plant and soil. The monthly maximal CO2 release and uptake rates range from 1.46 μmol·m−2 s−1 in October to 4.54 μmol·m−2 s−1 in April and from −16.62 μmol·m−2 s−1 in July to −2.59 μmol·m−2 s−1 in October. Diurnal peak Fc in this study is smaller than −22.0 μmol·m−2 s−1 for rice fields in Philippines [61] and −18.5 μmol·m−2 s−1 for tropical lowland flooded rice ecosystem in India [62].The seasonal and interannual variation in the monthly average Fc values from 2007 to 2018 are shown in Figure8(e) and the annual cumulative net CO2 ecosystem exchange (NEE) for each year is shown in Figure 10. Similar to Hs, the Fc has a distinct seasonal change during 2007–2013, demonstrating that it is higher during the transitional period (between April and May and between October and November) than that at other times and usually reaches the minimum value in July. Note that there is no obvious seasonal variation in Fc during 2014–2018. The annual NEE shows distinct jumps in magnitude, rapidly increasing after 2014 and varying widely from −966.9 to −75.6 g·C·m−2·yr−1, with a 12-year average value and a standard deviation of −522.0 and 295.7 g·C·m−2, respectively. This is possibly because of the changes of the underlying surface. The vegetation surrounding the PBL flux tower was altered from broad bean and rice to mixing plants including trees, and there was no crop rotation after 2014. Different plants also contribute to a large change in CO2 release and uptake. It is suggested that recording the changes of the underlying surface are necessary for the long-term eddy covariance measurements.Figure 10
Cumulative net CO2 ecosystem exchange (NEE) for each year during 2007 to 2018.


#### 3.3.2. Controlling Variables of CO2 Flux at Different Temporal Scales

The correlation coefficients of the Fc and environmental variables from half-hourly to yearly timescales are also given in Table6. Significantly negative correlation is found between Fc and Rn (r varying from −0.58 to −0.45) from half-hourly to monthly timescales. Similar results have been reported by Wang et al. [63] for a maize cropland in China. They also noticed that the daily Rn has a strong direct impact on Fc. During the day, CO2 absorption from plant photosynthesis increases with increasing Rn, which results in Fc becoming more negative [62]. Ta is a significantly negative relationship with Fc (r ranging from −0.41 to −0.27) between half-hourly and monthly timescales. The monthly average Fc is negatively correlated with the PPT and NDVI (r = −0.44 and −0.52, respectively). Like Hs and LE, WS is most significantly related to the total Fc on a yearly timescale (r = −0.92), which may be because of the changes of vegetation types and cover as well as the climate effect of local circulation. Overall, environmental drivers have varying influence on CO2 exchange at different timescales. This is coincident with several studies reported earlier [32, 49, 64].Table 6
Similar to Table4, the flux variable is the CO2 flux (Fc), and the environmental variables include Rn, Ta, WS, PPT, NDVI, and soil water content (SWC). The Fc is the average value from half-hourly to monthly timescales, and the annual Fc is the annual total value.
TimescalesEnvironmental variablesRnTaWSSWCPPTNDVIHalf-hourly−0.58∗−0.27∗−0.06∗−0.13∗Daily−0.45∗−0.33∗0.18∗−0.35∗−0.05∗Monthly−0.50∗−0.41∗0.12−0.45∗−0.44∗−0.52∗Yearly−0.40−0.20−0.92∗−0.380.060.13


## 3.1. Meteorological Conditions and Vegetation Cover

The climate at the measurement site features a monsoon climate of the low latitude plateau. Based on the observational data at an automatic weather station (AWS) located approximately 220 m southwest of the study site, the monthly average Ta during 1951–2018 ranged from 8.6°C in January to 20.2°C in June, with an annual average Ta of 15.1°C (Figure4). The monthly average PPT during 1951–2018 ranged from 12.9 mm in December to 218.4 mm in August, with the annual average PPT of 1048.1 mm and nearly 85% of precipitation falling in May to October. The monthly average WS during 1951–2018 ranged from 1.4 m·s−1 in September to 3.7 m·s−1 in March with the annual average WS of 2.4 m·s−1 and the strongest WS of 40.8 m·s−1. The gale is frequently experienced from winter through spring.Figure 4
Box chart of monthly averaged air temperature (a) and precipitation (b) from 1951 to 2018. The circle, square, and five-pointed star represent the mean, minimum, and maximum values, respectively. The horizontal line on the box represents the median and the box represents the 25–75th percentiles. The whisker represents the 5–95th percentile.
(a)(b)The meteorological variables during the study period are given in Table3 and shown in Figure 5. The difference of every variable in the dry and wet season is significant, namely, Rg, Ta, and RH, in the wet season are higher than those in the dry season. The solar energy resources are rich, with the annual integrated Rg varying from 5066 to 6703 MJ·m−2·yr−1, and a 12-year average value of 6154.1 MJ·m−2. The daily average Ta at heights of 2, 4, 10, and 20 m are always above 0°C, with the maximum and minimum values recorded in June 2014 and December 2013. The daily average Ts at the upper soil layer are higher than those at the lower soil layer in the wet season; conversely, Ts at the lower soil layer are higher than those at the upper soil layer in the dry season. The monthly average Ts values reach the maximum in July and the lowest in January. The RH within the near-surface layer usually decreases with increasing measurement height. The monthly average RH values reach the highest in August and the lowest in February. The WS within the near-surface layer commonly increases with increasing measurement height. The monthly average WS values reach the strongest in February and the weakest in September. The SWC values respond significantly to PPT. With the beginning of the wet season, the SWC rapidly increases and reaches its maximum.Table 3
Annual total global solar radiation (Rg), annual averaged wind speed (WS) at 10 m height, air temperature (Ta), and relative humidity (RH) at 4 m height, vapor pressure deficit (VPD), soil temperature (Ts), and soil water content (SWC) at 4 cm depth, and annual total precipitation (PPT) at the Dali site from 2007 to 2018.
YearRg (MJ·m−2)WS (m·s−1)Ta (°C)RH (%)VPD (kPa)Ts (°C)SWC (m−3·m−3)PPT (mm)20075066.12.616.668.00.6618.00.3201139.820085659.32.515.868.30.6218.40.3811322.820096702.72.716.462.10.7617.30.3521031.220106043.62.816.362.60.7416.40.394863.820116126.72.515.762.40.7215.90.426740.120126120.72.716.457.70.8315.80.422853.220136154.72.415.957.60.8416.10.423732.920146687.12.416.558.00.8515.90.332812.920156440.72.016.161.00.8115.90.326984.920166069.81.915.964.80.7015.30.3741153.620176356.92.016.163.60.7316.40.381909.820186421.61.915.563.00.7216.90.3441000.7Average6154.12.416.162.40.7516.50.373962.1Standard deviation450.20.30.33.60.070.90.038179.3Figure 5
Daily integrated global solar radiation (Rg) (a), daily averaged air temperature (Ta) (b), relative humidity (RH) (c), and wind speed (WS) (d) at four heights, daily averaged soil temperature (Ts) at five depths (e), and daily averaged soil water content (SWC) at five depths and daily total precipitation (PPT) (f) from 2007 to 2018.
(a)(b)(c)(d)(e)(f)The maximum daily total PPT value is observed in October 2015 (121.0 mm·d−1), and the maximum monthly total PPT is observed in June 2008 (336.7 mm·month−1). The annual total PPT fluctuates from 732.9 mm·yr−1 (2013) to 1322.8 mm·yr−1 (2008), with an average value of 962.1 mm (Table 3). The PPT in the wet season ranges from 643.2 to 1171.9 mm, accounting for 75.9% to 93.6% of the annual total PPT. Because of the frequent drought influence, the PPT for five consecutive years from 2010 to 2014 is below the 12-year average value.Local complicated terrain often affects the airflow within the atmospheric boundary layer, which not only changes wind speed but also changes wind direction. The predominant wind direction is easterly and east-southeasterly throughout the year (FigureS1). The prevailing and subprevailing winds at 10 m height display a clear diurnal cycle, with east-southeasterly and easterly winds during the daytime and west-northwesterly and static winds at night (Figure S2). There is a good correlation between the time when the prevailing winds switches and the sunrise and sunset (solar radiation heating effect). This result may be chiefly caused by the orientation of the CM and Erhai Lake, which extends from northwest to southeast.The NDVI has an obvious change before and after November 2013 because of the surface alternation (Figure6). For each year, the NDVI has two peaks in January or February (the fastest growing period of broad bean) and in July or August (the fastest growing period of rice), and it has also two valleys in April or May and in September or October (the harvesting periods of broad bean and rice) from 2007 to 2013. However, the NDVI has only a peak in July or August (growing period of the mixing plants) and a valley in January or February (nongrowing period of the mixing plants) from 2014 to 2018.Figure 6
Seasonal and interannual variations of the normalized difference vegetation index (NDVI) at 250 m spatial resolution and 16-day temporal resolution from 2007 to 2018.


## 3.2. Surface Energy Fluxes

### 3.2.1. Diurnal, Seasonal, and Interannual Variations in Surface Energy Fluxes

The monthly average diurnal variations in the energy balance components from 2007 to 2018 are shown in Figures7(a)–7(d). As expected, both heat and evapotranspiration vary with available net radiation. For each month, the Rn, Hs, LE, and Gs all exhibit clear, similar diurnal courses. During the day, Rn mainly comes from solar shortwave radiation and remains positive values, reaching its peak when the sun is in the zenith (between 12 : 00 and 13 : 00, depending on the season), while Rn chiefly comes from surface longwave radiation and becomes a negative value at night, reaching its minimum between 19 : 00 and 21 : 00. The monthly valleys and peaks of hourly average Rn vary from −67 W·m−2 in February to −32 W·m−2 in October and from 418 W·m−2 in December to 513 W·m−2 in June, respectively. Hs also attains its maximum and minimum when Rn reaches the peak and valley, which is consistent with the results found by Acharya et al. [47] for an agricultural field in Nepal. The monthly maximal and minimal Hs values vary from 66 W·m−2 in July to 127 W·m−2 in February and from −52 W·m−2 in March to −7 W·m−2 in August. The peak times of both LE and Gs are approximately 1 hour later than those of Rn and Hs. Similarly, Ajao et al. [48] reported that there was a lag of about 1 hour between Hs and LE at an agricultural site in Nigeria due to a phase change (water to vapor). LE remains a positive value and implies evaporation throughout the day. The monthly maximal and minimal LE values vary from 205 W·m−2 in December to 337 W·m−2 in August and from 3 W·m−2 in January to 20 W·m−2 in July. The diurnal average Gs is relatively small throughout the year, with its valley occurring in November (−25 W·m−2) and the peak occurring in June (47 W·m−2). The diurnal variations of Hs and LE in this study is obviously different from that observed in Erhai Lake [49], which may be related to the difference between land and water surface.Figure 7
Monthly average diurnal variations in net radiation flux (Rn) (a), sensible heat flux (Hs) (b), latent heat flux (LE) (c), and soil heat flux (Gs) (d), and monthly average diurnal variation in CO2 flux (Fc) (e) from 2007 to 2018. The error bars indicate plus or minus one standard deviation.
(a)(b)(c)(d)(e)The seasonal and interannual variations in the monthly average energy balance components from 2007 to 2018 are shown in Figures8(a)–8(d) and 9, respectively. For each year, the Rn, LE, and Gs values in the wet season are all larger than those in the dry season. The monthly average Rn generally attains the peak in June or July and the valley in January or December, with the maximum and minimum values of 229.0 and 22.8 W·m−2. The monthly average Hs shows a significant seasonal variation from 2007 to 2013, while Hs remains high and does not present evident variability from 2014 to 2018. Meanwhile, the annual average Hs increases from approximately 6 W·m−2 during 2007–2012 to 19 W·m−2 during 2013–2018. Although LE shows strong changes between dry and wet seasons during the entire period, the seasonal change amplitude of LE becomes small after 2014, and the annual average LE decreases from approximately 110 W·m−2 during 2007–2013 to 79 W·m−2 during 2014–2018. Overall, there is an obvious step change for Hs after 2013 and for LE after 2014, which may mainly be regulated by the vegetation types and cover around the PBL flux tower.Figure 8
Seasonal and interannual variations in monthly averaged net radiation flux (Rn) (a), sensible heat flux (Hs) (b), latent heat flux (LE) (c), soil heat flux (Gs) (d), and CO2 flux (Fc) (e) from 2007 to 2018.
(a)(b)(c)(d)(e)Figure 9
Annual averaged sensible heat flux (Hs), latent heat flux (LE), net radiation flux (Rn), and soil heat flux (Gs) from 2007 to 2018.


### 3.2.2. Surface Energy Partition

Owing to the change of vegetation cover, the proportions of the sensible and latent heat fluxes in net radiation (Hs/Rn and LE/Rn) are distinctly different. The annual average Hs/Rn increased from 5.3% during 2007–2012 to 17.7% during 2013–2018, while the annual average LE/Rn decreased from 94.6% during 2007–2013 to 72.8% during 2014–2018. In addition, the annual average Hs/Rn and LE/Rn during 2007–2018 are 11.2% and 85.8%, respectively. The annual average Bowen ratio (Hs/LE) during the study period is 0.15, suggesting that the surface heat exchange is modulated by the LE in the whole year. These findings are consistent with the study of Wang et al. [50] who reported that LE was the main part of energy consumption, with its annual value being more than twice of Hs in the karst region of the YGP, but not consistent with Li et al. [7] who found that the major energy source providing heat to the atmosphere was the Hs during the premonsoon period but the LE during the monsoon period at four flux stations (Linzhi, Nagqu, NamCo, and Qomolangma) over the TP. Through this comparative analysis, it is implied that the characteristics of energy exchange over the TP, YGP, and their transitional zones are similar in the wet season but clearly different in the dry season due to the influence of monsoon [7], the difference of location, and land-cover type.Furthermore, the monthly average Gs usually reaches the peak in May or June and the valley in November or December (Figure7(d)). The annual average proportion of the soil heat flux in net radiation (Gs/Rn) decreased from −0.09% during 2007–2013 to −6.3% during 2014–2018.


### 3.2.3. Controlling Variables in Hs at Different Temporal Scales

The energy, H2O and CO2 exchanges between the atmosphere and croplands depend not only on biological control but also on environmental factors such as meteorological conditions, vegetation characteristics, soil properties, and management practices [51–54]. Relationships between the Hs and environmental variables on timescales from half-hourly to yearly are given in Table 4. Compared to other environmental variables, the WS is the most important controlling variable for Hs from daily to yearly timescales. The correlation coefficients between them become higher as the timescales become longer, with r ranging from −0.95 to −0.55. The temperature difference between the upper soil layer (at 4 cm depth) and the atmosphere (ΔT) as well as WS multiplied by ΔT present positive effects on Hs on daily and monthly timescales but negative effects on Hs at half-hourly and yearly timescales. The correlation coefficients between Hs and WS are higher than those between Hs and ΔT as well as WS multiplied by ΔT, especially on the yearly timescale. Rn has significantly positive effects on Hs from half-hourly to monthly timescales (r ranging from 0.26 to 0.83) but a negative effect at yearly timescale (r = −0.39). The correlation coefficients between Hs and the rest of the environmental factors are relatively low.Table 4
Correlation coefficients of the averaged sensible heat flux (Hs) and environmental variables including the averaged wind speed (WS) at 4 m height, temperature difference between the upper soil layer (at 4 cm depth) and the atmosphere (ΔT), the WS multiplied by ΔT (WS∗ ΔT), net radiation (Rn) at 1.5 m height, air temperature (Ta), the total precipitation (PPT), and the normalized difference vegetation index (NDVI) at multiple timescales from 2007 to 2018. There is no examination about the PPT at half-hourly timescales and the NDVI between half-hourly and daily timescales.
TimescalesEnvironmental variablesWSΔTWS∗ ΔTRnTaPPTNDVIHalf-hourly0.01∗−0.20∗−0.10∗0.83∗0.31∗Daily−0.55∗0.19∗0.28∗0.53∗0.17∗−0.13∗Monthly−0.75∗0.38∗0.48∗0.26∗0.27∗0.25∗−0.08Yearly−0.95∗−0.39−0.32−0.39−0.49−0.100.03Note: the asterisk (∗) indicates the significance level of 0.05.There may be several reasons for the very high correlation of Hs with WS at the annual timescale. On the one hand, the alteration of the vegetation cover results in the surface roughness change, which is high enough to affect significantly the WS, that is, the average WS of 2.0 m·s−1 during 2014 to 2018 is significantly lower than that of 2.6 m·s−1 during 2007 to 2013 (Table 3). At the same time, the Hs significantly increases (Figure 8(b)). The close correlation between WS and Hs was also reported in many other studies. On the other hand, it is well known that the Hs is primarily determined by ΔT as well as by the turbulent exchange coefficient [55]. There is a strong local circulation within the boundary layer existing in the study area, which is caused by the thermally generated airflow (lake-land breeze and mountain-valley breeze), and the dynamically generated airflow (gorge wind) [56, 57]. This local circulation could modify wind speed and direction [33] and result in the air column to be vertically well mixed [58]. It may create more turbulent eddies and promote the vertical exchange of heat, H2O, and CO2 between the atmosphere and cropland when the WS increases but is less than a threshold. Wind-driven mixing largely modifies meteorological factors, which results in the changes of heat and water vapor exchange [59].


### 3.2.4. Controlling Variables in LE at Different Temporal Scales

The correlation coefficients of the LE and its environmental controls at multiple timescales are also summarized in Table5. Relative to other environmental variables, the product of WS and VPD is the main environmental variable controlling LE, showing a significantly positive correlation over all timescales, with r ranging from 0.41 to 0.72. Although WS and VPD are positively correlated with LE from half-hourly to monthly timescales (r being more than 0.32), the values of r between WS and LE are slightly smaller than those between VPD and LE. On a yearly timescale, WS has a strong positive correlation with LE (r = 0.95), while VPD has a weak negative correlation with LE (r = −0.10). This implies that the effect of WS is not significant as VPD on LE when the timescale is below the monthly timescale, whereas WS plays a decisive role in affecting LE on a yearly scale. Clearly, the impacts of WS and VPD, and their products on LE relies on timescales [60]. There are positive correlations between Rn and LE on all timescales with r ranging from 0.31 to 0.82, which indicates that Rn primarily regulates the variation in water vapor flux, particularly on half-hourly timescales. Ta displays significant positive correlations with LE below monthly timescales but nonsignificant correlation at the yearly timescale.Table 5
Similar to Table4, the flux variable is the latent heat flux (LE), and the environmental variables include the WS, Rn, Ta, PPT, NDVI, as well as the vapor pressure deficit (VPD), and the WS multiplied by VPD (WS ∗ VPD) at 4 m height.
TimescalesEnvironmental variablesWSVPDWS∗VPDRnTaPPTNDVIHalf-hourly0.32∗0.53∗0.43∗0.82∗0.49∗Daily0.32∗0.45∗0.41∗0.66∗0.56∗−0.11∗Monthly0.42∗0.43∗0.45∗0.58∗0.61∗0.21∗0.11Yearly0.95∗−0.100.72∗0.310.34−0.08−0.12


## 3.2.1. Diurnal, Seasonal, and Interannual Variations in Surface Energy Fluxes

The monthly average diurnal variations in the energy balance components from 2007 to 2018 are shown in Figures7(a)–7(d). As expected, both heat and evapotranspiration vary with available net radiation. For each month, the Rn, Hs, LE, and Gs all exhibit clear, similar diurnal courses. During the day, Rn mainly comes from solar shortwave radiation and remains positive values, reaching its peak when the sun is in the zenith (between 12 : 00 and 13 : 00, depending on the season), while Rn chiefly comes from surface longwave radiation and becomes a negative value at night, reaching its minimum between 19 : 00 and 21 : 00. The monthly valleys and peaks of hourly average Rn vary from −67 W·m−2 in February to −32 W·m−2 in October and from 418 W·m−2 in December to 513 W·m−2 in June, respectively. Hs also attains its maximum and minimum when Rn reaches the peak and valley, which is consistent with the results found by Acharya et al. [47] for an agricultural field in Nepal. The monthly maximal and minimal Hs values vary from 66 W·m−2 in July to 127 W·m−2 in February and from −52 W·m−2 in March to −7 W·m−2 in August. The peak times of both LE and Gs are approximately 1 hour later than those of Rn and Hs. Similarly, Ajao et al. [48] reported that there was a lag of about 1 hour between Hs and LE at an agricultural site in Nigeria due to a phase change (water to vapor). LE remains a positive value and implies evaporation throughout the day. The monthly maximal and minimal LE values vary from 205 W·m−2 in December to 337 W·m−2 in August and from 3 W·m−2 in January to 20 W·m−2 in July. The diurnal average Gs is relatively small throughout the year, with its valley occurring in November (−25 W·m−2) and the peak occurring in June (47 W·m−2). The diurnal variations of Hs and LE in this study is obviously different from that observed in Erhai Lake [49], which may be related to the difference between land and water surface.Figure 7
Monthly average diurnal variations in net radiation flux (Rn) (a), sensible heat flux (Hs) (b), latent heat flux (LE) (c), and soil heat flux (Gs) (d), and monthly average diurnal variation in CO2 flux (Fc) (e) from 2007 to 2018. The error bars indicate plus or minus one standard deviation.
(a)(b)(c)(d)(e)The seasonal and interannual variations in the monthly average energy balance components from 2007 to 2018 are shown in Figures8(a)–8(d) and 9, respectively. For each year, the Rn, LE, and Gs values in the wet season are all larger than those in the dry season. The monthly average Rn generally attains the peak in June or July and the valley in January or December, with the maximum and minimum values of 229.0 and 22.8 W·m−2. The monthly average Hs shows a significant seasonal variation from 2007 to 2013, while Hs remains high and does not present evident variability from 2014 to 2018. Meanwhile, the annual average Hs increases from approximately 6 W·m−2 during 2007–2012 to 19 W·m−2 during 2013–2018. Although LE shows strong changes between dry and wet seasons during the entire period, the seasonal change amplitude of LE becomes small after 2014, and the annual average LE decreases from approximately 110 W·m−2 during 2007–2013 to 79 W·m−2 during 2014–2018. Overall, there is an obvious step change for Hs after 2013 and for LE after 2014, which may mainly be regulated by the vegetation types and cover around the PBL flux tower.Figure 8
Seasonal and interannual variations in monthly averaged net radiation flux (Rn) (a), sensible heat flux (Hs) (b), latent heat flux (LE) (c), soil heat flux (Gs) (d), and CO2 flux (Fc) (e) from 2007 to 2018.
(a)(b)(c)(d)(e)Figure 9
Annual averaged sensible heat flux (Hs), latent heat flux (LE), net radiation flux (Rn), and soil heat flux (Gs) from 2007 to 2018.


## 3.2.2. Surface Energy Partition

Owing to the change of vegetation cover, the proportions of the sensible and latent heat fluxes in net radiation (Hs/Rn and LE/Rn) are distinctly different. The annual average Hs/Rn increased from 5.3% during 2007–2012 to 17.7% during 2013–2018, while the annual average LE/Rn decreased from 94.6% during 2007–2013 to 72.8% during 2014–2018. In addition, the annual average Hs/Rn and LE/Rn during 2007–2018 are 11.2% and 85.8%, respectively. The annual average Bowen ratio (Hs/LE) during the study period is 0.15, suggesting that the surface heat exchange is modulated by the LE in the whole year. These findings are consistent with the study of Wang et al. [50] who reported that LE was the main part of energy consumption, with its annual value being more than twice of Hs in the karst region of the YGP, but not consistent with Li et al. [7] who found that the major energy source providing heat to the atmosphere was the Hs during the premonsoon period but the LE during the monsoon period at four flux stations (Linzhi, Nagqu, NamCo, and Qomolangma) over the TP. Through this comparative analysis, it is implied that the characteristics of energy exchange over the TP, YGP, and their transitional zones are similar in the wet season but clearly different in the dry season due to the influence of monsoon [7], the difference of location, and land-cover type.Furthermore, the monthly average Gs usually reaches the peak in May or June and the valley in November or December (Figure7(d)). The annual average proportion of the soil heat flux in net radiation (Gs/Rn) decreased from −0.09% during 2007–2013 to −6.3% during 2014–2018.


## 3.2.3. Controlling Variables in Hs at Different Temporal Scales

The energy, H2O and CO2 exchanges between the atmosphere and croplands depend not only on biological control but also on environmental factors such as meteorological conditions, vegetation characteristics, soil properties, and management practices [51–54]. Relationships between the Hs and environmental variables on timescales from half-hourly to yearly are given in Table 4. Compared to other environmental variables, the WS is the most important controlling variable for Hs from daily to yearly timescales. The correlation coefficients between them become higher as the timescales become longer, with r ranging from −0.95 to −0.55. The temperature difference between the upper soil layer (at 4 cm depth) and the atmosphere (ΔT) as well as WS multiplied by ΔT present positive effects on Hs on daily and monthly timescales but negative effects on Hs at half-hourly and yearly timescales. The correlation coefficients between Hs and WS are higher than those between Hs and ΔT as well as WS multiplied by ΔT, especially on the yearly timescale. Rn has significantly positive effects on Hs from half-hourly to monthly timescales (r ranging from 0.26 to 0.83) but a negative effect at yearly timescale (r = −0.39). The correlation coefficients between Hs and the rest of the environmental factors are relatively low.Table 4
Correlation coefficients of the averaged sensible heat flux (Hs) and environmental variables including the averaged wind speed (WS) at 4 m height, temperature difference between the upper soil layer (at 4 cm depth) and the atmosphere (ΔT), the WS multiplied by ΔT (WS∗ ΔT), net radiation (Rn) at 1.5 m height, air temperature (Ta), the total precipitation (PPT), and the normalized difference vegetation index (NDVI) at multiple timescales from 2007 to 2018. There is no examination about the PPT at half-hourly timescales and the NDVI between half-hourly and daily timescales.
TimescalesEnvironmental variablesWSΔTWS∗ ΔTRnTaPPTNDVIHalf-hourly0.01∗−0.20∗−0.10∗0.83∗0.31∗Daily−0.55∗0.19∗0.28∗0.53∗0.17∗−0.13∗Monthly−0.75∗0.38∗0.48∗0.26∗0.27∗0.25∗−0.08Yearly−0.95∗−0.39−0.32−0.39−0.49−0.100.03Note: the asterisk (∗) indicates the significance level of 0.05.There may be several reasons for the very high correlation of Hs with WS at the annual timescale. On the one hand, the alteration of the vegetation cover results in the surface roughness change, which is high enough to affect significantly the WS, that is, the average WS of 2.0 m·s−1 during 2014 to 2018 is significantly lower than that of 2.6 m·s−1 during 2007 to 2013 (Table 3). At the same time, the Hs significantly increases (Figure 8(b)). The close correlation between WS and Hs was also reported in many other studies. On the other hand, it is well known that the Hs is primarily determined by ΔT as well as by the turbulent exchange coefficient [55]. There is a strong local circulation within the boundary layer existing in the study area, which is caused by the thermally generated airflow (lake-land breeze and mountain-valley breeze), and the dynamically generated airflow (gorge wind) [56, 57]. This local circulation could modify wind speed and direction [33] and result in the air column to be vertically well mixed [58]. It may create more turbulent eddies and promote the vertical exchange of heat, H2O, and CO2 between the atmosphere and cropland when the WS increases but is less than a threshold. Wind-driven mixing largely modifies meteorological factors, which results in the changes of heat and water vapor exchange [59].


## 3.2.4. Controlling Variables in LE at Different Temporal Scales

The correlation coefficients of the LE and its environmental controls at multiple timescales are also summarized in Table5. Relative to other environmental variables, the product of WS and VPD is the main environmental variable controlling LE, showing a significantly positive correlation over all timescales, with r ranging from 0.41 to 0.72. Although WS and VPD are positively correlated with LE from half-hourly to monthly timescales (r being more than 0.32), the values of r between WS and LE are slightly smaller than those between VPD and LE. On a yearly timescale, WS has a strong positive correlation with LE (r = 0.95), while VPD has a weak negative correlation with LE (r = −0.10). This implies that the effect of WS is not significant as VPD on LE when the timescale is below the monthly timescale, whereas WS plays a decisive role in affecting LE on a yearly scale. Clearly, the impacts of WS and VPD, and their products on LE relies on timescales [60]. There are positive correlations between Rn and LE on all timescales with r ranging from 0.31 to 0.82, which indicates that Rn primarily regulates the variation in water vapor flux, particularly on half-hourly timescales. Ta displays significant positive correlations with LE below monthly timescales but nonsignificant correlation at the yearly timescale.Table 5
Similar to Table4, the flux variable is the latent heat flux (LE), and the environmental variables include the WS, Rn, Ta, PPT, NDVI, as well as the vapor pressure deficit (VPD), and the WS multiplied by VPD (WS ∗ VPD) at 4 m height.
TimescalesEnvironmental variablesWSVPDWS∗VPDRnTaPPTNDVIHalf-hourly0.32∗0.53∗0.43∗0.82∗0.49∗Daily0.32∗0.45∗0.41∗0.66∗0.56∗−0.11∗Monthly0.42∗0.43∗0.45∗0.58∗0.61∗0.21∗0.11Yearly0.95∗−0.100.72∗0.310.34−0.08−0.12


## 3.3. CO2 Flux

### 3.3.1. Diurnal, Seasonal, and Interannual Variations in CO2 Flux

The monthly average diurnal variations in Fc from 2007 to 2018 are shown in Figure7(e). For each month, the Fc has a noticeable diurnal cycle with positive values at night reaching their maximum at approximately 07 : 00 and negative values during the day reaching their minimum at approximately 13 : 00. This result indicates that the study area acts as a weak carbon source at night because of ecosystem respiration, but acts as a carbon sink during the day as the CO2 uptake sourced from plant photosynthesis significantly exceeds CO2 emission released from plant and soil. The monthly maximal CO2 release and uptake rates range from 1.46 μmol·m−2 s−1 in October to 4.54 μmol·m−2 s−1 in April and from −16.62 μmol·m−2 s−1 in July to −2.59 μmol·m−2 s−1 in October. Diurnal peak Fc in this study is smaller than −22.0 μmol·m−2 s−1 for rice fields in Philippines [61] and −18.5 μmol·m−2 s−1 for tropical lowland flooded rice ecosystem in India [62].The seasonal and interannual variation in the monthly average Fc values from 2007 to 2018 are shown in Figure8(e) and the annual cumulative net CO2 ecosystem exchange (NEE) for each year is shown in Figure 10. Similar to Hs, the Fc has a distinct seasonal change during 2007–2013, demonstrating that it is higher during the transitional period (between April and May and between October and November) than that at other times and usually reaches the minimum value in July. Note that there is no obvious seasonal variation in Fc during 2014–2018. The annual NEE shows distinct jumps in magnitude, rapidly increasing after 2014 and varying widely from −966.9 to −75.6 g·C·m−2·yr−1, with a 12-year average value and a standard deviation of −522.0 and 295.7 g·C·m−2, respectively. This is possibly because of the changes of the underlying surface. The vegetation surrounding the PBL flux tower was altered from broad bean and rice to mixing plants including trees, and there was no crop rotation after 2014. Different plants also contribute to a large change in CO2 release and uptake. It is suggested that recording the changes of the underlying surface are necessary for the long-term eddy covariance measurements.Figure 10
Cumulative net CO2 ecosystem exchange (NEE) for each year during 2007 to 2018.


### 3.3.2. Controlling Variables of CO2 Flux at Different Temporal Scales

The correlation coefficients of the Fc and environmental variables from half-hourly to yearly timescales are also given in Table6. Significantly negative correlation is found between Fc and Rn (r varying from −0.58 to −0.45) from half-hourly to monthly timescales. Similar results have been reported by Wang et al. [63] for a maize cropland in China. They also noticed that the daily Rn has a strong direct impact on Fc. During the day, CO2 absorption from plant photosynthesis increases with increasing Rn, which results in Fc becoming more negative [62]. Ta is a significantly negative relationship with Fc (r ranging from −0.41 to −0.27) between half-hourly and monthly timescales. The monthly average Fc is negatively correlated with the PPT and NDVI (r = −0.44 and −0.52, respectively). Like Hs and LE, WS is most significantly related to the total Fc on a yearly timescale (r = −0.92), which may be because of the changes of vegetation types and cover as well as the climate effect of local circulation. Overall, environmental drivers have varying influence on CO2 exchange at different timescales. This is coincident with several studies reported earlier [32, 49, 64].Table 6
Similar to Table4, the flux variable is the CO2 flux (Fc), and the environmental variables include Rn, Ta, WS, PPT, NDVI, and soil water content (SWC). The Fc is the average value from half-hourly to monthly timescales, and the annual Fc is the annual total value.
TimescalesEnvironmental variablesRnTaWSSWCPPTNDVIHalf-hourly−0.58∗−0.27∗−0.06∗−0.13∗Daily−0.45∗−0.33∗0.18∗−0.35∗−0.05∗Monthly−0.50∗−0.41∗0.12−0.45∗−0.44∗−0.52∗Yearly−0.40−0.20−0.92∗−0.380.060.13


## 3.3.1. Diurnal, Seasonal, and Interannual Variations in CO2 Flux

The monthly average diurnal variations in Fc from 2007 to 2018 are shown in Figure7(e). For each month, the Fc has a noticeable diurnal cycle with positive values at night reaching their maximum at approximately 07 : 00 and negative values during the day reaching their minimum at approximately 13 : 00. This result indicates that the study area acts as a weak carbon source at night because of ecosystem respiration, but acts as a carbon sink during the day as the CO2 uptake sourced from plant photosynthesis significantly exceeds CO2 emission released from plant and soil. The monthly maximal CO2 release and uptake rates range from 1.46 μmol·m−2 s−1 in October to 4.54 μmol·m−2 s−1 in April and from −16.62 μmol·m−2 s−1 in July to −2.59 μmol·m−2 s−1 in October. Diurnal peak Fc in this study is smaller than −22.0 μmol·m−2 s−1 for rice fields in Philippines [61] and −18.5 μmol·m−2 s−1 for tropical lowland flooded rice ecosystem in India [62].The seasonal and interannual variation in the monthly average Fc values from 2007 to 2018 are shown in Figure8(e) and the annual cumulative net CO2 ecosystem exchange (NEE) for each year is shown in Figure 10. Similar to Hs, the Fc has a distinct seasonal change during 2007–2013, demonstrating that it is higher during the transitional period (between April and May and between October and November) than that at other times and usually reaches the minimum value in July. Note that there is no obvious seasonal variation in Fc during 2014–2018. The annual NEE shows distinct jumps in magnitude, rapidly increasing after 2014 and varying widely from −966.9 to −75.6 g·C·m−2·yr−1, with a 12-year average value and a standard deviation of −522.0 and 295.7 g·C·m−2, respectively. This is possibly because of the changes of the underlying surface. The vegetation surrounding the PBL flux tower was altered from broad bean and rice to mixing plants including trees, and there was no crop rotation after 2014. Different plants also contribute to a large change in CO2 release and uptake. It is suggested that recording the changes of the underlying surface are necessary for the long-term eddy covariance measurements.Figure 10
Cumulative net CO2 ecosystem exchange (NEE) for each year during 2007 to 2018.


## 3.3.2. Controlling Variables of CO2 Flux at Different Temporal Scales

The correlation coefficients of the Fc and environmental variables from half-hourly to yearly timescales are also given in Table6. Significantly negative correlation is found between Fc and Rn (r varying from −0.58 to −0.45) from half-hourly to monthly timescales. Similar results have been reported by Wang et al. [63] for a maize cropland in China. They also noticed that the daily Rn has a strong direct impact on Fc. During the day, CO2 absorption from plant photosynthesis increases with increasing Rn, which results in Fc becoming more negative [62]. Ta is a significantly negative relationship with Fc (r ranging from −0.41 to −0.27) between half-hourly and monthly timescales. The monthly average Fc is negatively correlated with the PPT and NDVI (r = −0.44 and −0.52, respectively). Like Hs and LE, WS is most significantly related to the total Fc on a yearly timescale (r = −0.92), which may be because of the changes of vegetation types and cover as well as the climate effect of local circulation. Overall, environmental drivers have varying influence on CO2 exchange at different timescales. This is coincident with several studies reported earlier [32, 49, 64].Table 6
Similar to Table4, the flux variable is the CO2 flux (Fc), and the environmental variables include Rn, Ta, WS, PPT, NDVI, and soil water content (SWC). The Fc is the average value from half-hourly to monthly timescales, and the annual Fc is the annual total value.
TimescalesEnvironmental variablesRnTaWSSWCPPTNDVIHalf-hourly−0.58∗−0.27∗−0.06∗−0.13∗Daily−0.45∗−0.33∗0.18∗−0.35∗−0.05∗Monthly−0.50∗−0.41∗0.12−0.45∗−0.44∗−0.52∗Yearly−0.40−0.20−0.92∗−0.380.060.13


## 4. Conclusions

We report the diurnal, seasonal, and interannual variability in energy and CO2 fluxes and their environmental drivers at various temporal scales from 2007 to 2018 at the Dali observational site over the southeast extension of the TP. The Rn, Hs, LE, and Gs all have similar diurnal courses, reaching their maximum values around noon and attaining their minimum values around early evening. Moreover, these factors present obvious seasonal changes, with their values during the wet season being larger than those during the dry season. The average values and standard deviations of Rn, Hs, LE, and Gs for 12 years are 112.8 ± 12.9, 12.7 ± 7.1, 96.8 ± 16.4, and −0.4 ± 0.5 W·m−2, respectively. The Fc has a noticeable diurnal cycle with positive values at night and negative values during the daytime, and it also exhibits clear seasonal changes with the maximum negative value in July and the highest positive value during the transitional periods between the dry and wet seasons. Due to the vegetation types and cover surrounding the PBL flux tower being altered after 2014, the annual total Fc values vary widely, fluctuating from −966.9 to −75.6 g·C·m−2·yr−1 with a 12-year average value and standard deviation of −522.0 ± 295.7 g·C·m−2. This result suggests that the study area acts as an overall carbon sink during all years.Meteorological variables have different effects on Hs, LE, and Fc at different temporal scales. The WS shows a decreasing trend according to the correlation coefficients with Hs above daily timescales, while the product of WS and VPD is more responsible for the variation of LE from the half-hourly to yearly timescales. For Fc, the Rn has high contributions at the daily to monthly scales, but the WS is the most important factor at the yearly scale. The turbulent exchanges of energy, H2O, and CO2 fluxes are largely influenced by the WS, particularly on the yearly scale. This phenomenon is caused by the alternation of vegetation cover surrounding the study site and is also likely caused by the climate effect of local circulation in the complex terrain region. The effects of WS on both energy and carbon exchange processes indicates the effect of topography in complex regions. More effort should be made about the atmosphere-land interaction, which will help to improve the numerical simulation in these regions.



---

*Source: 4329199-2022-12-23.xml*