## Abstract

The transfer process and variation characteristics of water and heat fluxes play an important role in both the ecosystem function and the climatic patterns of arid and semi-arid areas. However, both water (latent heat flux, Fe) and heat (sensible heat flux, Fh) fluxes, and the factors impacting each, operate on different time scales. Therefore, the wavelet transformation was used to reconstruct water and heat flux data measured by an eddy covariance system installed in Horqin sandy land, China. The distribution and response regularity of water and heat fluxes in response to control factors were then analysed across different time scales. The results indicated that the responses of Fe and Fh to controlling factors were similar on 30 min and daily scales, to which the energy terms of net radiation (Fn) and available energy (Fa) showed an extremely significant positive contribution. Removing the interference of the energy terms, Fe and Fh are strongly correlated with the vapour pressure deficit (VPD) and air temperature (Ta) at the daily scale. The water conditions of VPD and relative humidity (RH) were the main influencing factors at the monthly scale. At the seasonal scale, Fe was positively correlated with Fa and Fn but negatively correlated with soil heat flux (Fg), RH and wind speed (WS). On the contrary, Fh was positively correlated with RH but negatively correlated with Fn. At the annual scale, Fe and Fh showed a significant positive correlation with all factors. By the increasing time scale, the cross-correlations between various hydroclimatic processes become better. Generally, at the daily and annual scales, there is a high correlation between Fe, Fh and the control factors, and a pronounced periodicity is exhibited. Another factor that is known to cause such effects is the so-called long-term persistence (or else long-range dependence or the Hurst phenomenon). These results suggest that the time scales of flux and impacting factors should be considered when undertaking water and heat flux analyses of sandy land and other ecosystems.

## HIGHLIGHTS

Latent and sensible heat fluxes in sandy land vary in response to various influencing factors at different time scales.

The same period latent heat flux was less than sensible heat flux, indicating that the energy was mainly used to heat the surface air and make the surface temperature rise rapidly.

Wavelet coherence of water heat flux and energy component showed highly consistent daily and annual common periods.

### Graphical Abstract

## INTRODUCTION

The formation and development of sandy dunes are the result of a process of degradation and even collapse of grassland ecosystems, which are the product of global ecosystem degradation (Ru *et al.* 2018; Wang *et al.* 2019). Grassland ecosystems in arid and semi-arid areas are vulnerable, and under the interacting effects of climate change and human activities, grassland desertification has become a worldwide ecological problem (Pnuma 1994; Wang *et al.* 2015; Li *et al.* 2017). Grassland degradation and desertification can change the carbon cycle of terrestrial ecosystems under the climatic condition of persistent or aggravated drought (Schimel *et al.* 2001; Reynolds *et al.* 2017). The carbon sink of savanna succession ecosystems produced by recovery from grassland desertification and vegetation degradation often decreases and may even be converted back into a carbon source (Kolb *et al.* 2013; Ferlan & Eler 2016). Studies have shown that overgrazing, unreasonable land reclamation and climate change are the main causes of grassland desertification (Liang *et al.* 2009). The water and heat fluxes in sandy land are closely related to the desertification process, but their variation characteristics and the response to meteorological factors are still unclear.

In terms of material and energy exchange, different ecosystems (forests, farmlands, grasslands or wetlands, for example) are generally regarded as relatively independent, discrete units (Piao *et al.* 2009), responding to different eco-hydrological and meteorological factors. Hydrological processes of arid and semi-arid areas are extremely sensitive to regional meteorological conditions (Schnorbus *et al.* 2014; Stagl *et al.* 2014). Therefore, the changes in meteorological factors will have direct or indirect, temporal and spatial impacts on the eco-hydrological processes (Liu *et al.* 2015; Huang *et al.* 2019a, 2019b). The dynamic process that hydrological elements change continuously or periodically over time has attracted widespread attention.

The key factors causing the vulnerability of ecological environment in sandy land are lower and uneven distribution of precipitation between seasons and years (Foley *et al.* 2003). Since the 1990s, the climate has shown a warm and humid trend, which has led to the large-scale rapid development of desertification (Huang *et al.* 2016). Desertification is a process of species diversity attenuation and an increase of hydro-meteorological disturbance (Fang *et al.* 2009). Therefore, the study of water exchange and energy balance between land and air in sandy ecosystems provides an important basis for determining the eco-hydrological effects upon the natural vegetation pattern in arid and semi-arid areas.

The relationship between water and heat restricts the growth and survival of vegetation of sandy land, and this fragile ecosystem is very sensitive to climate change (Bai *et al.* 2004). Under the influence of climate change and human activities, the original water and heat cycle process of sand dunes have changed, which further affect the regional climate change and destroy the already fragile sand dune ecosystem (Law *et al.* 2002). The most obvious anthropogenic signal in the hydrological cycle appears to be the over-exploitation of groundwater, which has a visible effect on accelerating the pace of drought and desertification (Koutsoyiannis 2019). Water and heat flux are important determinants of both water vapour and energy exchange between land surface and atmosphere, which together constitute the driving mechanism of atmospheric circulation (Wang *et al.* 2016).

The exchange of water and heat flux in sandy land ecosystems is controlled by energy conditions, moisture conditions and aerodynamic conditions (Timouk *et al.* 2009; Duveiller *et al.* 2018; Che *et al.* 2019). The changes of these control conditions across multiple time scales lead to the different distribution characteristics of water and heat flux (Katul & Parlange 1995; Gruber *et al.* 2002). Due to the effects of multiple factors on water and heat balance, the temporal variability of mass and energy fluxes (Baldocchi *et al.* 2001) and the direct effect of land surface processes upon the climate can occur on time scales ranging from seconds to the seasonal and annual (Viterbo 2003). Therefore, it is of great significance to study the temporal variability of water and heat flux in relation to their controlling factors across different time scales in sandy land of arid and semi-arid areas.

To analyse the variation characteristics of water and heat fluxes over different time scales within a sand dune ecosystem, and to reveal its relationship with meteorological factors, this paper first shows the variation characteristics of meteorological elements and energy components over time. Then, the water and heat flux data observed by the eddy covariance (EC) system are transformed and reconstructed. Further, the characteristics of water and heat fluxes in sandy dunes are analysed in different time scales. Finally, the co-relationships between meteorological elements with water and heat fluxes across different time scales were discussed, and the primary factors influencing water and heat fluxes across different time scales were also identified.

## MATERIALS AND METHODS

### The study site

^{2}of Horqin sandy land (HSL) (118 °35′ to 123 °30′E, 42 °41′ to 45 °15′N), which is located in the Mongolian Plateau of north China (Figure 1). It was grassland in history and began to devolve into sandy land over 100 years ago (Han

*et al.*2010; Li

*et al.*2017). The experiment site is located within the HSL, and the area is 55 km

^{2}(122 °33′00″E to 122 °41′00″E, 43 °18′48″N to 43 °21′24″N) (Duan

*et al.*2011). It has a monsoon-influenced, temperate semi-arid continental climate, with an average annual precipitation of 389 mm, of which 69.3% falls during the growing season (June to August). The annual potential evapotranspiration (PET) is 1,412 mm (Tong

*et al.*2019). It has an average annual temperature of around 6.6 °C, and an average annual wind speed (WS) of 3.6 m·s

^{−1}. The prevalent wind direction is northwest in winter and spring, and southwest during summer and fall (Duan

*et al.*2011; Tong

*et al.*2019). The main types of sandy dunes in the study area are fixed dunes and semi-fixed dunes. Zonal and non-zonal soils are distributed alternately (Tong

*et al.*2016).

Since 2012, the experimental station has been set up on sandy dune lands, and the observation system was composed of micrometeorology, radiation and eddy correlation systems. The station was equipped with an open-path EC observation system, wind, temperature and humidity gradient sensors, which were installed at different heights (Table 1). Additional temperature and humidity sensors were used to correct the air temperature and humidity measured by the EC system. The height of EC system is 4.61 m and the sampling frequency is 10 Hz.

Meteorological elements . | Instruments and sensor . | Location . |
---|---|---|

CO_{2} and H_{2}O concentrations | LI-7500, LI-COR, USA | 4.61 m |

Vertical wind speed | CSAT3, CSI, USA | 4.61 m |

Net radiation | CNR-LITE, CSI, USA | 4.61 m |

Wind speed and direction | 014A and 034A-L, CSI, USA | 2 m, 10 m |

Air temperature and humidity | HMP45C, CSI, US | 0 m, 2 m |

Soil temperature | 105T, CSI, USA | 10 cm, 20 cm, 40 cm, 80 cm |

Soil moisture content | TDR, CS615, CSI, USA | 120 cm, 160 cm |

Soil heat flux | HFP01, CSI, USA | 10 cm, 20 cm, 30 cm |

Meteorological elements . | Instruments and sensor . | Location . |
---|---|---|

CO_{2} and H_{2}O concentrations | LI-7500, LI-COR, USA | 4.61 m |

Vertical wind speed | CSAT3, CSI, USA | 4.61 m |

Net radiation | CNR-LITE, CSI, USA | 4.61 m |

Wind speed and direction | 014A and 034A-L, CSI, USA | 2 m, 10 m |

Air temperature and humidity | HMP45C, CSI, US | 0 m, 2 m |

Soil temperature | 105T, CSI, USA | 10 cm, 20 cm, 40 cm, 80 cm |

Soil moisture content | TDR, CS615, CSI, USA | 120 cm, 160 cm |

Soil heat flux | HFP01, CSI, USA | 10 cm, 20 cm, 30 cm |

### EC data processing

In this paper, the format conversion (TOB3 to TOB1) of the 10 Hz raw data downloaded directly from the OPEC was carried out, and then a series of data corrections, such as outlier elimination and tilt correction, were carried out by using EddyPro, version 6.0 (LI-COR Biosciences, USA) (Falge *et al.* 2001) to obtain the flux data with a half-hour period. The half-hour data obtained were used for quality control such as the stationarity test and the turbulence characteristic test, and more accurate flux data (Fh, Fe, CO_{2}, H_{2}O) were obtained. Then the missing and deleted data were interpolated by the Australian Flux Network OzFluxQC system (Isaac *et al.* 2016). Through the above process, the continuous water and heat flux data without gaps were produced. Finally, the complete series of half-hour meteorological and flux data from 2013 to 2016 were used to carry out the wavelet transform (WT) analysis to explore the variation characteristics of water heat flux across different time scales.

We have provided more descriptive statistical information of site data to better reflect the statistical characteristics of time-series data. Table 2 shows the first four statistical moments (i.e., mean, standard deviation, skewness, and kurtosis coefficients) for every process (radiation flux, temperature, wind, humidity, etc.).

Statistics . | Fn . | Fg . | Fa . | Fe . | Fh . | Ta . | Ts . | VPD . | WS . | RH . |
---|---|---|---|---|---|---|---|---|---|---|

W·m^{−2}
. | W·m^{−2}
. | W·m^{−2}
. | W·m^{−2}
. | W·m^{−2}
. | °C . | °C . | kPa . | m·s^{−1}
. | % . | |

Mean | 51.696 | 0.906 | 51.363 | 19.927 | 19.718 | 7.347 | 9.712 | −0.278 | 3.594 | 66.518 |

SD | 150.003 | 30.015 | 130.885 | 35.902 | 58.367 | 14.029 | 13.479 | 1.106 | 2.162 | 23.714 |

Skewness | 1.340 | 1.444 | 1.267 | 3.744 | 1.418 | −0.153 | −0.028 | −0.102 | 0.758 | −0.190 |

Kurtosis | 0.775 | 2.633 | 0.655 | 0.051 | 3.490 | −1.219 | −1.246 | 1.215 | 0.442 | −1.036 |

Statistics . | Fn . | Fg . | Fa . | Fe . | Fh . | Ta . | Ts . | VPD . | WS . | RH . |
---|---|---|---|---|---|---|---|---|---|---|

W·m^{−2}
. | W·m^{−2}
. | W·m^{−2}
. | W·m^{−2}
. | W·m^{−2}
. | °C . | °C . | kPa . | m·s^{−1}
. | % . | |

Mean | 51.696 | 0.906 | 51.363 | 19.927 | 19.718 | 7.347 | 9.712 | −0.278 | 3.594 | 66.518 |

SD | 150.003 | 30.015 | 130.885 | 35.902 | 58.367 | 14.029 | 13.479 | 1.106 | 2.162 | 23.714 |

Skewness | 1.340 | 1.444 | 1.267 | 3.744 | 1.418 | −0.153 | −0.028 | −0.102 | 0.758 | −0.190 |

Kurtosis | 0.775 | 2.633 | 0.655 | 0.051 | 3.490 | −1.219 | −1.246 | 1.215 | 0.442 | −1.036 |

### The Mann–Kendall trend test of hydro-meteorological processes

*X*= {

*x*

_{1},

*x*

_{2}, …,

*x*}, the test statistic is given by:where and are the ranks of observations and of the time series, respectively. As can be seen from Equation (2), the test statistic depends only on the ranks of the observations, rather than their actual values, resulting in a distribution-free test statistic. Distribution-free tests have the advantage that their power and significance are not affected by the actual distribution of the data.

_{n}*S*statistic in Equation (1) above are given by:where

*n*is the number of observations.

*S*tends to normality as the number of observations becomes large. The significance of trends can be tested by comparing the standardized variable

*Z*in Equation (5) with the standard normal variate at the desired significance level

*a*, where the subtraction or addition of unity in Equation (5) is a continuity correction (Kendall 1975).

We mainly focus on *Z* and *P* values. If the *Z* value is less than 0, it indicates a downward trend, while if the *Z* value is greater than 0, it indicates an upward trend. In the process of our study (Table 3), *Z* (Fa, Fe, Fh, Ta, Ts, and VPD) is greater than 0, *P* < 0.001; *Z* (Fn) is greater than 0, *P* < 0.05; *Z* (WS) is greater than 0, *P* < 0.01; *Z* (Fg, RH) is less than 0, *P* < 0.001. Thus, we reject the null hypothesis (assuming that there is no trend), indicating that there is a statistically significant trend change in the sequence data.

### Wavelet transform and correlation analysis

*et al.*2008). The WT is a kind of signal processing technology, which is used to quantify the spectral characteristics of data that vary with time and space (Baldocchi & Wilson 2001; Baldocchi

*et al.*2001), and has been applied more and more in the analysis of unstable time series. The WT is defined as:where the WT has two variables: scale () and translation ().

*a*controls the stretch out and draw back, and controls the translation of wavelet function. The scale corresponds to the frequency (inverse ratio), and the amount of translation corresponds to time. The complex conjugate of the scaled and translated mother wavelet is given as .

In this paper, the continuous wavelet transform (CWT) and the wavelet transform coherence (WTC) toolboxes for MATLAB are used to reflect the correlation between two time series. The CWT can be viewed as a microscope capable of revealing the details of the signal at the position with the magnification *a*. This approach is useful for time-series analysis, where the power of the signal is modulated at different scales which vary with time (Vargas *et al.* 2010). CWT assumes that the data are circular, while for the actual time series a non-stationary model may found more appropriate. Therefore, due to the influence of edge effect, there are some errors at the beginning and the end of wavelet power spectrum (WPS). The cone of influence (COI) is defined as the area affected by the edge (Torrence & Compo 1998). In the area outside the COI, the influence of the edge effect can be ignored.

*et al.*2004), and the multi-scale (hourly, daily, monthly, seasonal, and annual) feature extraction is carried out from the flux data with a period of half-hourly. The Morlet wavelet is given by the following equation:where (the base frequency) is set to 6. The Morlet is a complex nonorthogonal wavelet with a good time and scale resolution. It has the advantage of returning information for both amplitude and phase.

The WTC is defined as a smoothing in time and scale. The smoothing is performed in both time and scale using the method proposed by Torrence & Webster (2000), which provides the minimal amount of smoothing necessary to include two independent points in both dimensions. WTC can be thought of as the local correlation between two time series. Importantly, it finds regions in time–frequency space where two time series covary but do not necessarily have a high common power (Grinsted *et al.* 2004; Cazelles *et al.* 2008). In order to compare the variation characteristics of the water and heat fluxes of different time series, the problem of missing data of the fluxes and the non-stationarity of the time series were accounted for (Stoy *et al.* 2005). Then, the WT was used to reconstruct the data of water and heat flux to different time scales, followed by the extraction of distribution characteristics for different time scales. Wavelet analysis can effectively solve the non-stationarity of time series, and it is an appropriate tool to explain the contribution of environmental variables to the time series of water and heat flux (Katul & Parlange 1995; Liu *et al.* 2017).

The calculation and testing of the Pearson correlation coefficient are based on the cor function and cor.test function of the corrplot package in R (Legendre & Gauthier 2014). The correlation analysis matrix is also obtained by the corrplot package using R software. Correlation analysis is used to analyse the relationship between two variables. In real analysis, correlation analysis often has the influence of the third or even more variables, so that the correlation coefficient cannot truly reflect its linear correlation degree. A partial correlation analysis was used to solve this problem. Partial correlation analysis refers to the process of eliminating the influence of these conditional variables when two variables are related to one or more other variables (conditional variables) at the same time and only analysing the correlation degree between the two variables. The judgment index is the *R* value of the correlation coefficient (Yuan *et al.* 2015).

### The climacogram and standardized cross-climacogram

Stochastic processes are families of random variables (denoted as , where underlined symbols denote random variables and *t* denotes time) that are often used to represent the temporal evolution of natural processes. Natural processes as well as their mathematical representation as stochastic processes evolve in continuous time (Dimitriadis & Koutsoyiannis 2015). The hydrological and meteorological processes in the real world change in all time and scales. In addition, these processes are interdependent. Dependence and change are closely related. Now, we may take a note that dependence should not be interpreted as memory, but as a change. In particular, long-term dependence is not a long-term memory but long-term change (Koutsoyiannis 2021). How is change quantified in stochastics? A better choice is to make this quantification in a stochastic way. In this case, we see change as variability across different time scales. Furthermore, variability is quantified according to variance.

Referring to the half-hourly records of 4-year (2013–2016) length time series extracted from stations of the EC system, which we denote as , , … , , we take the following steps:

- 1.
- 2.
- 3.
We repeat the same procedure to form time series at time scales 4, 6, 8, … , up to scale 17,280 (one-fourth of the record length) and calculate the variance , , , … , .

- 4.
We plot (in double logarithmic axes) the variance as a function of time scale

*k*.

The function of the variance vs. time scale is called the climacogram (Koutsoyiannis 2010). According to the above steps, then we have the empirical climacogram . If we have assumed a model for these processes and we determine the variance, , from the model, we have the theoretical climacogram. Notably, if we have produced the times series from a model, the empirical climacogram will not necessarily coincide with the theoretical, because there is estimation bias. To make them coincide, we must subtract the bias from the theoretical climacogram (Koutsoyiannis 2021).

*H*− 2, where

*H*is the so-called Hurst parameter which takes on values in the interval (0, 1). The case where this slope is constant for all time scales corresponds to a simple scaling behaviour (e.g. Koutsoyiannis 2006):which defines the Hurst–Kolmogorov (HK) process (Kolmogorov 1940; Hurst 1951).

The superscript (*D*) in denotes the time step of discretization.

*k*of the averaged process, we can define the cross-climacogram (Koutsoyiannis 2019):where and

*h*is lag.

*k*and lag

*h*to further expand to describe the dependence of different processes (take the dependence between water and heat flux (Fe and Fh) and energy component (Fn) as an example).

## RESULTS AND DISCUSSION

### Meteorological and energy components characteristics

^{−1}. Ta, Ts and VPD all showed an obvious single peak in each year. The peak of Ta and Ts appeared in July of each year, while the peak of VPD appeared in May, June, July and August, across years, respectively. However, this trend was not obvious in WS, which showed a state of continual fluctuation (Figure 2(a) and 2(b)).

The energy component is shown in Figure 2(c), with an obvious seasonal pattern and an interannual variation. The change of Fn, Fa (Fn–Fg) and Fg fluctuated between 0 to 500 W·m^{−2}, −50 to 70 W·m^{−2}, and 0 to 400 W·m^{−2}, respectively. Fn and Fa showed a single peak every year, and the peak appeared in May (2013), July (2014), July (2015), and June (2016). The characteristic of Fg was stronger in summer and weaker in winter.

The value of Fe was less than Fh in the same period, which indicated that the energy was mainly used to heat the surface air and make the surface temperature rise rapidly in the study area. Fh decreased successively in 2013, 2015, 2014, and 2016, while Fe decreased year by year in 2013, 2015, 2016, and 2014. In short, except for 2014, the changes in each year were relatively consistent, which showed a single peak curve on the whole. Fe in 2014 showed a significant difference from other years. There was no obvious trend of rising first and then falling in the whole year, but always showed a fluctuating state, and there was only a high state in August (Figure 2(d)).

The diurnal variation of water and heat fluxes in the sand dune ecosystem was significant, with a single peak for both Fh and Fe. This is consistent with the results of previous studies (Wang *et al.* 2009, 2016; Huang *et al.* 2019a, 2019b), and this change is mainly affected by solar height, sky cloud cover and micro-meteorological conditions nearer to the earth's surface. Time series of the water and heat fluxes exhibits pronounced annual cycles, which are primarily explained by the modulation of the yearly flux amplitude by other variables, such as the Fn and Fa (Brunsell & Cassandra 2013).

We apply this climacogram method on the multiple time-scale indices of half-hourly time series. Figure 1 depicts the climacograms of meteorological, energy components and water heat flux of sandy land. In these cases, the observed behaviour is spectacularly different from white noise, while the Hurst–Kolmogorov behaviour is evident with the Hurst parameter *H* as high as 0.98 for Ta, Ts, and VPD. The slope for large scales again suggests a strongly persistent behaviour with the Hurst parameter *H* = 0.82–0.98 for above hydroclimatic processes.

*H*= 1/2 corresponds to white noise as the slope is −1. High values of

*H*(>1/2) indicate an enhanced change at large scales, else known as long-term persistence, or strong clustering (grouping) of similar values. This is quite common in hydroclimatic processes. Low values

*H*(<1/2) indicate quite a different behaviour, called antipersistence. Such behaviour is much less frequent in hydroclimatic processes. In Figure 3, the last point of the climacogram falls outside the power-law fit; however, we ignore these last scales, since this is due to the estimation of variance from a small sample at these scales.

### Variation of water and heat fluxes in different time scales

*et al.*2020).

By comparing and analysing the correlation coefficients of Fe, Fh and each control factor in different time scales, it was found that the correlation is generally higher at the daily and annual scales, which is consistent with the variations in diurnal and annual periods shown by the wavelet coherence spectrum (Figure 4), which may be caused by the rotation and revolution of the earth (Xu *et al.* 2005).

In general, the local correlation was higher at Fh than at Fe, suggesting a strong link between variations in the water heat flux and the energy component (Fn, Fa), especially in 2014. There were four clusters on the annual periodicity, and each cluster gathers in June, July, and August. But in 2014, there was only a short period of aggregation of Fe in August, and there was no obvious fluctuation state in other time scales (between time 2.3 and 3.5 × 10^{4} min in Figure 4(1) and 4(2)). The cause for emergence is that the proportion of Fe in Fn in 2014 was much smaller than that of Fh. This may have been caused by the existence of the flux tower in the study area, which makes the sandy vegetation canopy more sparse, causing an uneven and relatively less seasonal distribution of precipitation in this year (Wilson *et al.* 2002).

### Response of water and heat fluxes to control factors

The clustering relationship visible in Figure 5 shows that the relationship between the elements is divided into three parts. Fa, Fn, Fg, Ta, and Ts belong to the energy part, where Ta and Ts belong to the temperature group. VPD belongs to the part of the water condition and WS belongs to the part of the aerodynamic condition. Further cluster analysis shows that Fa, Fn, Fg, Fh and Fe belong to the same group.

. | Fn . | Fg . | Fa . | Fe . | Fh . | Ta . | Ts . | VPD . | WS . | RH . |
---|---|---|---|---|---|---|---|---|---|---|

Z | 1.607 | −36.45 | 16.009 | 16.332 | 17.445 | 14.675 | 25.033 | 23.878 | 2.7581 | −5.087 |

P | 0.0107 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00581 | 0.000000364 |

Trend | 0.44 | −0.32 | 0.60 | 1.21 | 2.90 | 0.17 | 0.32 | 0.05 | 0.02 | −0.87 |

. | Fn . | Fg . | Fa . | Fe . | Fh . | Ta . | Ts . | VPD . | WS . | RH . |
---|---|---|---|---|---|---|---|---|---|---|

Z | 1.607 | −36.45 | 16.009 | 16.332 | 17.445 | 14.675 | 25.033 | 23.878 | 2.7581 | −5.087 |

P | 0.0107 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00000000000000022 | 0.00581 | 0.000000364 |

Trend | 0.44 | −0.32 | 0.60 | 1.21 | 2.90 | 0.17 | 0.32 | 0.05 | 0.02 | −0.87 |

The Pearson correlation coefficients (*r*) of different time scales are shown in Table 4.

. | Fe . | Fh . | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

30 min . | Hourly . | Daily . | Monthly . | Seasonal . | Annual . | 30 min . | Hourly . | Daily . | Monthly . | Seasonal . | Annual . | |

Fa | 0.50** | 0.05 | 0.73** | 0.18* | 0.62** | 0.98*** | 0.80*** | − 0.05 | 0.91*** | 0.44** | − 0.16* | 0.99*** |

Fg | 0.36** | − 0.08 | 0.68** | − 0.19* | − 0.33** | 0.80*** | 0.65** | 0.08 | 0.80*** | 0.15* | 0.08 | 0.88*** |

Fn | 0.51** | 0.03 | 0.76** | 0.19* | 0.52** | 0.96*** | 0.81*** | 0.03 | 0.93*** | 0.35** | − 0.27* | 0.99*** |

RH | − 0.03 | − 0.04 | − 0.46** | 0.30** | − 0.39** | 0.83*** | − 0.14* | 0.01 | − 0.46** | − 0.46** | 0.25* | 0.74** |

Ta | 0.38** | 0.01 | 0.53** | − 0.14* | 0.16* | 0.98*** | 0.37** | 0.02 | 0.69** | 0.09 | − 0.08 | 0.95*** |

Ts | 0.36** | − 0.05 | 0.46** | − 0.05 | 0.34** | 0.96*** | 0.33** | 0.06 | 0.51** | 0.18* | − 0.15* | 0.91*** |

VPD | 0.36** | 0.07 | 0.57** | − 0.25* | 0.04 | 0.91*** | 0.35** | 0.04 | 0.67** | 0.19* | − 0.12* | 0.88*** |

WS | 0.16* | 0.04 | 0.61** | 0.14* | − 0.37** | 0.03 | 0.26* | 0.03 | 0.75** | 0.21* | 0.17* | 0.18 |

. | Fe . | Fh . | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

30 min . | Hourly . | Daily . | Monthly . | Seasonal . | Annual . | 30 min . | Hourly . | Daily . | Monthly . | Seasonal . | Annual . | |

Fa | 0.50** | 0.05 | 0.73** | 0.18* | 0.62** | 0.98*** | 0.80*** | − 0.05 | 0.91*** | 0.44** | − 0.16* | 0.99*** |

Fg | 0.36** | − 0.08 | 0.68** | − 0.19* | − 0.33** | 0.80*** | 0.65** | 0.08 | 0.80*** | 0.15* | 0.08 | 0.88*** |

Fn | 0.51** | 0.03 | 0.76** | 0.19* | 0.52** | 0.96*** | 0.81*** | 0.03 | 0.93*** | 0.35** | − 0.27* | 0.99*** |

RH | − 0.03 | − 0.04 | − 0.46** | 0.30** | − 0.39** | 0.83*** | − 0.14* | 0.01 | − 0.46** | − 0.46** | 0.25* | 0.74** |

Ta | 0.38** | 0.01 | 0.53** | − 0.14* | 0.16* | 0.98*** | 0.37** | 0.02 | 0.69** | 0.09 | − 0.08 | 0.95*** |

Ts | 0.36** | − 0.05 | 0.46** | − 0.05 | 0.34** | 0.96*** | 0.33** | 0.06 | 0.51** | 0.18* | − 0.15* | 0.91*** |

VPD | 0.36** | 0.07 | 0.57** | − 0.25* | 0.04 | 0.91*** | 0.35** | 0.04 | 0.67** | 0.19* | − 0.12* | 0.88*** |

WS | 0.16* | 0.04 | 0.61** | 0.14* | − 0.37** | 0.03 | 0.26* | 0.03 | 0.75** | 0.21* | 0.17* | 0.18 |

*Note*: Pearson test significant level *P-*value (0.001, 0.01, 0.05), corresponding symbols (‘***’, ‘**’, ‘*’).

The correlations between Fe and its controlling factors were significantly varied at different time scales. At an hourly scale, none of the correlation coefficients were significant. However, at the daily scale, the correlation between Fe and all control factors was significant (*P* < 0.01). Among them, the highest value of *r* was 0.76 (Fe and Fn), while Fn was mainly related to Fa, so the value of *r* between Fe and Fa was also high. Fe and RH showed a negative correlation, with an *r* value of −0.46. At the monthly scale, the highest *r* occurred between Fe and RH at 0.3 (*P* < 0.01), while the lowest *r* was appeared between Fe and VPD at −0.25 (*P* < 0.05). At the seasonal scale, the highest *r* occurred between Fe and Fa (0.62, *P* < 0.01) followed by Fe and Fn (0.52, *P* < 0.01). The *r* was then negative between Fe and both RH and WS (−0.39 and −0.37, *P* < 0.01). At the annual scale, there was a very significant positive correlation between Fe and all control factors (*P* < 0.001). The highest *r* appeared between Fe and both Fa and Ta, both were 0.98 followed by Fn and Ts, which was 0.96, while WS was the only factor not significant with Fe.

Fe- and energy-related factors such as Fa, Fn, Ta, and Ts showed significant correlations at 30 min, daily, seasonal and annual scales but lower correlations in hourly and monthly scales. However, Fe and Fg showed a negative correlation at the monthly and seasonal scales, especially at the seasonal scale where it possessed a significant negative correlation. Fe was significantly positively correlated at the annual scale with both the water-related factors RH and VPD, but significantly negatively correlated with RH and positively correlated with VPD at the daily scale. In addition, Fe also showed positive and negative correlations with RH at the monthly and seasonal scales, respectively, but was not significantly correlated with VPD. The significant correlation between Fe and the aerodynamic factor, WS, was significant mainly at the daily and seasonal scales, but there was a positive correlation at the daily scale and a negative correlation at the seasonal scale.

Similar to Fe, Fh also was differently correlated with each control factor at different time scales. At the hourly scale, there was no significant correlation between Fh and each control factor. At the daily scale, *r* between Fh and Fn was the highest (0.93), while Fn was mainly related to Fa, so *r* between Fh and Fa was also relatively high (0.91, *P* < 0.001). Fh and RH showed a negative correlation, and *r* was −0.46 (*P* < 0.01). At the monthly scale, the highest *r* occurred between Fh and Fa which was 0.44 (*P* < 0.01), while the lowest *r* was appeared between Fh and RH at −0.46 (*P* < 0.01). However, the seasonal scale was opposite to the monthly scale with the highest *r* occurring between Fh and RH (0.25, *P* < 0.05) and a negative *r* between Fh and Fn (−0.27, *P* < 0.05). At the annual scale, Fh showed a very significant positive correlation with each control factor (*P* < 0.001). The highest *r* value of Fh with Fa and Fn was 0.99 for both, followed by Ta and Ts, which were 0.95 and 0.91, respectively.

In summary, Fh showed significant positive correlations with the energy-related factors (Fa, Fg, Fn, Ta, and Ts) at the 30 min, daily and annual scales, but lower correlations at the hourly, monthly, and seasonal scales, and mostly negative correlations at the seasonal scale. A linear or segmental linear correlation was shown between Fe, Fh, and the energy-component factors. Solar radiation was the main source of surface energy and the basis of matter and energy exchange (Liang *et al.* 2010), so energy factors play a leading role in the effects of Fe and Fh. For the water-related factors RH and VPD, Fh and RH showed significant positive correlation at the annual scale but significant negative correlation at the daily and monthly scales, while Fh showed a significant positive correlation with VPD at the daily and annual scales but a negative correlation at the seasonal scale. At the monthly scale, the response of Fe and Fh to RH differed completely. Overall, RH affects Fe and Fh to a lesser degree, with conflicting trends. At the seasonal scale, the correlation between Fe and each control factor was generally higher than Fh. However, Fh showed an abnormal response compared to other time scales, exhibiting a positive correlation with RH and a negative correlation with Fn at the seasonal level. This suggests that Fh is negatively limited by Fn on this time scale. This may be because the seasonal allocation of Fn to Fh is strongly influenced by soil moisture content (Hao *et al.* 2007), which is inhibited when soil moisture reaches the threshold of stress. However, vegetation phenology may also be a major factor leading to flux partitioning at the seasonal scale in the study area (Jiang *et al.* 2015).

*ρ*≈ 0.5, while for positive lags

*ρ*> 0.5, Fe with a maximum

*ρ*= 0.55 at the scale of 30 min and at a lag of 12 h, Fh with a maximum

*ρ*= 0.55 at the scale of 30 min and at a lag of 24 h. This does not look too high because the cross-correlations between the latent and sensible heat fluxes with eight various hydroclimatic processes are weak.

We can reduce the process variance by averaging at larger time scales. Then the correlation becomes stronger, reaching a peak of *ρ* = 0.73 (Fe) and 0.68 (Fh) for time scale of 1 year, beyond which the maximum *ρ* decreases. Thus, the evidence for causality becomes stronger by increasing the time scale up to a certain point. By the increasing time scale, the cross-correlations between various hydroclimatic processes become better.

### The influencing factors of different time scales

This paper further investigates the differences in the response of Fh and Fe to each control factor at different time scales of hourly, daily, monthly, seasonal, and annual on the basis of half-hourly data (Figure 7). At the daily and annual scales, the water and heat fluxes were mainly influenced by meteorological elements and were significantly correlated with almost all meteorological factors. While at other scales, water and heat fluxes were influenced by a single factor, or by several factors combined. Here, there may be other relevant factors in addition to the meteorological ones, such as vegetation status, especially at monthly and seasonal scales (Wu *et al.* 2018a).

Because the water heat fluxes (Fe and Fh) and energy components belong to the same energy group, the relationship between them may mask the influence of environmental variables on water and heat fluxes. Therefore, this study seeks to further understand the mechanism of daily and annual environmental factors influencing water and heat fluxes through partial correlation analysis. We found that on the daily scale, excluding the interference of the energy component, there was a significant partial correlation between water heat flux and VPD (Table 5), which was stronger than for other environmental factors. Fe and Fh were negatively correlated with RH, but their inhibitory effect was extremely weak. VPD can directly stimulate vegetation to reduce evapotranspiration by closing stomatal openings, but this response may take a long time to be detected (Yu *et al.* 2013). At the same time, the diurnal variation of temperature makes the water and heat fluxes as well as the VPD of sandy ecosystems more periodic at the daily scale (Ouyang *et al.* 2014). The strong correlation between water heat flux and VPD on the daily scale may be caused by the joint influence of these two mechanisms. At daily time scales, interactions between physiological and biophysical features of the canopy and its micro-climates are the dominant mode of flux variation (Katul *et al.* 2001). RH, however, negatively affects both to a lesser extent, with a higher relative humidity leading to a slower energy exchange and a smaller the change in water heat flux.

Exclude variables . | Independent variables . | Response variables . | |||
---|---|---|---|---|---|

Fe . | Fh . | ||||

Daily . | Annual . | Daily . | Annual . | ||

Energy components (Fn, Fa, Fg) | RH | − 0.31* | 0.51** | − 0.27* | 0.97*** |

Ta | 0.39* | 0.73** | 0.44** | 0.79** | |

Ts | 0.12 | 0.64** | 0.21* | 0.73** | |

VPD | 0.52** | 0.63** | 0.67** | 0.61** | |

WS | 0.14 | − 0.84** | 0.18 | − 0.77** |

Exclude variables . | Independent variables . | Response variables . | |||
---|---|---|---|---|---|

Fe . | Fh . | ||||

Daily . | Annual . | Daily . | Annual . | ||

Energy components (Fn, Fa, Fg) | RH | − 0.31* | 0.51** | − 0.27* | 0.97*** |

Ta | 0.39* | 0.73** | 0.44** | 0.79** | |

Ts | 0.12 | 0.64** | 0.21* | 0.73** | |

VPD | 0.52** | 0.63** | 0.67** | 0.61** | |

WS | 0.14 | − 0.84** | 0.18 | − 0.77** |

*Note*: Pearson test significant level *P-*value (0.001, 0.01, 0.05), corresponding symbols (‘***’, ‘**’, ‘*’).

At the annual scale, Fe and Fh are significantly correlated with each environmental factor, among them, with a highly significant linear positive correlation with Ta and negative correlation with WS. Since Ta affects the magnitude of saturated water pressure in air and WS affects the speed of water vapour diffusion on the seasonal patterns (Xiao *et al.* 2018), Ta and WS also directly affect the variation of Fe and Fh. Particularly, the correlations of Fe and Fh with WS are extremely weak with the exception of annual time scales. In most situations, winds can accelerate evaporation and provide the dynamic conditions for convective exchange, but when the wind speed reaches a certain threshold, the extent of its effect on water heat flux will be reduced (Krishnan *et al.* 2012) even to the point of having a negative effect. Just as at the annual scale, Fe and WS are negatively correlated, i.e., wind speed has an inhibitory effect on Fe. In general, at longer time scales, seasonal patterns share importance with an interannual climate change (Katul *et al.* 2001).

Interestingly, the response characteristics of water heat fluxes to the control factors differed significantly at different time scales, especially at the monthly and seasonal scales. RH inhibited Fh at the monthly scale and facilitated it at the seasonal scale; conversely, RH facilitated Fe at the monthly scale and restricted it at the seasonal scale. On the monthly and seasonal scales, RH is likely to have an indirect impact on water heat flux by affecting evapotranspiration. The effect of RH on water heat flux may also be regulated by soil moisture, which may be lower in sandy land during our study period. This provides a consideration for future further research on water heat flux in different months and growing seasons.

## CONCLUSIONS

Using wavelet analysis, this paper investigates the characteristics of influencing factors at different time scales upon water heat fluxes in typical sandy ecosystems and draws the following conclusions.

The annual variations of energy components (Fn, Fa, and Fg) and meteorological elements (Ta, Ts, and VPD) were consistent, showing a single peak trend, which were stronger in summer and weaker in spring and winter. Over time, Fh was dominant relative to Fe of the same period in Fn, which indicated that Fn in this area is mainly used to heat the surface air, resulted in significant temperature difference between the surface and the atmosphere, thus accelerating the heat exchange between them through conduction and convection.

The response characteristics of Fh and Fe to the control factors in sandy land ecosystem at multiple time scales differed significantly, and there are still many uncertainties at the monthly and seasonal scales. At the daily and annual scales, the correlations between Fe, Fh, and the control factors were generally higher, and intensity varied in diurnal and annual periods. Therefore, we detail understandings of the characteristics of water and heat fluxes at different time scales to better explore the fluxes.

## ACKNOWLEDGEMENTS

This research was supported by the National Natural Science Foundation of China (Nos. 52079063 and 51779116) and the Natural Science Foundation of Inner Mongolia Autonomous Region of China (2019JQ06).

## DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

*&*Management