Time scale and influencing factor variations of water and heat fluxes in sandy land

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.

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.

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₂, H₂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.).

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.

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.

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).

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.
  We calculate the variance estimate ⁠, where ‘1’ indicates the time scale of half-hourly as:

  

  (8)

  where is the estimate of the mean and n = 70,127 is the sample size.
• 2.
  We form a time series at time scale 2 (half-hours) by averaging pairs of consecutive items of the time series, i.e.,

  

  (9)

  and we calculate the estimate of the variance in a similar manner.
• 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).

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

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⁻², −50 to 70 W·m⁻², and 0 to 400 W·m⁻², 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.

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⁴ 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).

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.

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

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).

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.

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.

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.

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.

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 cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.