Evapotranspiration dynamics and their drivers in a temperate mixed forest in northeast China

Evapotranspiration (ET) is a vital part of the global water cycle and is closely related to carbon sequestration. Analysing ET dynamics and their drivers would benefit for improving our understanding of the global water and carbon cycles. Using an eddy covariance (EC) approach, we analysed ET dynamics and their drivers in a temperate mixed forest over northeast China from 2016 to 2017. The results showed that 43.55% of our eddy covariance data passed the quality control. In addition, the energy balance ratio was 0.62, indicating that measurements were reliable. The measured ET showed clear single peak patterns with seasonal and diurnal variations. The daily ET ranged from 0 to 7.75 mm d−1 and the hourly ET ranged from 0 to 0.28 mm h−1. The ranges of hourly ET floated from 0 to 0.05 mm h−1 at non-growing season (November to April) while ranged from 0 to 0.28 mm h−1 at active growing season (May to October). The diurnal ET dynamics during the non-growing season were driven by air temperature (Ta), but were governed by global radiation (Rg) during the active growing season. Leaf area index (LAI) comprehensively reflected the variations of Ta and Rg, and was found to be the primary factor shaping the seasonal dynamics of ET. The annual ET rates were 501.91 ± 5.30 mm year−1 and 554.60 ± 11.24 mm year−1 for 2016 and 2017, respectively. Therefore, energy supply, represented by Ta and Rg, governed ET dynamics in our temperate mixed forest, while variables representing the energy supply affecting ET dynamics differed among seasons and time scales. ET dynamics indicated that a temperate mixed forest is important to the global water cycle. Our results improved our understanding of ET dynamics in the studied region.


INTRODUCTION
Evapotranspiration (ET) is defined as the water lost as vapour rising from the land to the atmosphere. It is a vital part of the global water cycle (Kool et al., 2014;Wang et al., 2010;Wang & Dickinson, 2012). ET is also closely associated with carbon sequestration as water loss and carbon sequestration both primarily appear in the stoma of a plant (Chapin, Matson & Mooney, 2012), which made ET closely associate with the global carbon cycle (Zhu et al., 2015b). Analysing the dynamics of ET and their drivers is beneficial for understanding the mechanisms underlying ET dynamics, thus accurately modeling ET. This also results in a better understanding of the global water and carbon cycles (Anapalli et al., 2020).
Analysing ET dynamics and their drivers requires accurate measurements (Li et al., 2009;Rana & Katerji, 2000;Wang & Dickinson, 2012). The micrometeorological theory known as eddy covariance (EC) has been widely used in measuring ET for its in-situ, continuous, and accurate measurements (Baldocchi & Ryu, 2011;Drexler et al., 2004;Tanaka et al., 2008). Many studies have analysed ET dynamics and their drivers for various ecosystems using EC (Cristiano et al., 2015;Kume et al., 2011;Ma et al., 2015;Song et al., 2017;Xu et al., 2014;Yue et al., 2019). Results showed that daily accumulated ET and their annual values differed among ecosystems (Kumagai et al., 2005;Li et al., 2010;Tong et al., 2017). Some studies have analysed ET dynamics and their drivers for temperate forests in North America (Schaefer et al., 2014;Wang et al., 2015;Wehr et al., 2017), Europe (Gielen et al., 2010;Pita et al., 2013;Soubie et al., 2016), and Asia (Tsuruta et al., 2016;Wu et al., 2013;Zhang et al., 2012). Carbon fluxes (Ohtsuka et al., 2005;Saigusa et al., 2002) and resource use efficiency (Okada, Takagi & Nishida, 2019) in temperate mixed forests over Asia have also been intensively reported, while less attention has been given to ET dynamics and what drives them in these ecosystems, which made the ET dynamics in this region poorly understood. In addition, temperate mixed forests over Asia are climax communities and the succession climax indicates the potential of this region to regulate regional energy and mass cycles (Zhang, Han & Yu, 2006;Zhang et al., 2006b). Analysing ET dynamics and their drivers in temperate mixed forests will improve our understanding of ET dynamics in this region, which benefited for revealing the potential role of this region in global water cycles.
We analysed ET dynamics and their drivers for a temperate mixed forest based on EC measurements in a broad-leaved Korean pine forest in northeast China. The main objectives of this study were: (1) to reveal ET dynamics for a temperate mixed forest in Northeast China, (2) to clarify the major drivers of ET dynamics,and (3) to reveal the magnitude of annual ET in this unique ecosystem. In order to address these objectives, we first addressed the energy balance closure of the eddy covariance measurements, which was directly related to the ET measurements (Wilson et al., 2002). We then showed the dynamics of ET and their environmental factors, including solar radiation (R g ), air temperature (T a ), relative humidity (RH), precipitation, and leaf area index (LAI) to determine the drivers of ET dynamics. Then the drivers of ET dynamics were performed. Our results elucidate the understanding of ET dynamics and ET modeling over a temperate mixed forest, which improves our understanding of the water cycle processes and their linked carbon cycles.

Site description
EC measurements were taken at the Yichun broad-leaved Korean pine forest experimental station (48.2292 N, 129.2661 E, 342 m.a.s.l) from January 2016 to December 2017. This station was located in the middle of Heilongjiang Province, which was the northeastern-most province of China (Fig. 1A). The experimental ecosystem had a flat topography with an evenly distributed mixed forest (14,141 hm −2 ). The mean fetch of the tower in all directions was lower than 400 m, with the prevailing southwest and northeast wind directions (Fig. 1B). The forest was primarily comprised of Korean Pine (Pinus koraiensis), Mongolian Oak (Quercus mongolica), Amur Linden (Tilia amurensis), and Maple Birch (Betula costata), with a density of 1,152 plants per hm −2 and an average stand age of 220-years (Fig. 1C). Korean Pine (Pinus koraiensis), Mongolian Oak (Quercus mongolica), Amur Linden (Tilia amurensis), and Maple Birch (Betula costata) accounted for 10-15% of all plants in the study area, respectively. The mean canopy height was 26 m with a maximum and minimum LAI of 6.2 and 0.5 m 2 m −2 , respectively. The experimental forest had a temperate continental climate with a long-term mean annual air temperature (MAT) of 0.63 C and a mean annual precipitation (MAP) of 610.7 mm. The highest T a appeared in July with a long-term mean  value of 20.4 C, while the lowest T a occurred in January with a value of −22.5 C. The greatest amount of precipitation occurred in July (147 mm), while the lowest precipitation occurred in February (6.7 mm). Precipitation was commonly seen as snowfall from November to April, and totaled approximately 90 mm year −1 . The soil in the experimental ecosystem was characterized as dark brown soil with high soil organic matter content.

Data measurements
We employed a flux tower with a height of 70 m to measure the ET and conventional meteorological variable dynamics.
The EC system was composed of a datalogger (Model CR5000; Campbell Scientific Inc., Logan, UT, USA), a 3-D sonic anemometer (Model CSAT3; Campbell Scientific Inc., Logan, UT, USA), and an infrared gas analyzer (Model LI-7500; Licor Inc., Lincoln, NE, USA). These were used to measure the net H 2 O exchange between the terrestrial ecosystem and the atmosphere, which can be considered as the ET. The EC sensors, including CSAT3 and LI-7500, were mounted at a height of 50 m, which was about twice of the canopy height. The EC measurement sampled the raw data with a frequency of 20 Hz and recorded the flux data at 30-min intervals.
Conventional meteorological variables were also measured, including air temperature (T a ), relative humidity (RH), soil temperature (T s ), precipitation, and soil heat flux (G). T a and RH were measured at different heights (1.5, 20, 30, 40, 50, 60 m) with shielded and aspirated probes (Model HMP45C; Campbell Scientific Inc., Logan, UT, USA). T s was measured using thermocouple probes (Model 105T; Campbell Scientific Inc., Logan, UT, USA) at different depths (0, 5, 10, 15, 20, 40, 80, 160, 320 cm). Precipitation was recorded with a rain gauge (Model 52203; R.M. Young Company, Traverse City, MI, USA) Radiation data, including global radiation (R g ), reflected radiation (R r ), downward long wave radiation (DLR), and upward long wave radiation (ULR), were measured with a radiation heat balance station (CAMS620-SP, Huatron, China) at 1-h intervals. The measured 1-h data were linearly interpolated into half-hour scale to match the time series of ET and other environmental variables. Net radiation (R n ) was calculated as the sum of net short-wave radiation (R ns ) and net long-wave radiation (R nl ), where R ns was the difference between R g and R r while R nl was the difference between DLR and ULR.

Flux data processing
The net H 2 O exchange, which was deemed as ET and latent heat (LE), was calculated as the covariance between vertical wind speed (v) fluctuations and water vapor concentration fluctuations. The spikes of wind speeds and H 2 O concentrations were detected and deleted based on CSAT3 measuring three-dimension wind speeds and LI-7500 measuring H 2 O concentrations (Vickers & Mahrt, 1997). The raw ET was deemed as missing if the missing or deleted data were more than 3,600 in each half-hour (10% of all measurements) (Wen et al., 2010). The raw ET was calculated as the covariance of vertical wind (v) speed fluctuations and H 2 O concentration fluctuations using the despiked data. The calculated raw fluxes were derived from a three-dimensional rotation (Aubinet et al., 2000), Webb, Pearman and Leuning (WPL) correction (Webb, Pearman & Leuning, 1980), spectral correction for high-frequency losses (Moncrieff et al., 1997), low-frequency losses (Moncrieff et al., 2005), physical instrument separation losses (Horst & Lenschow, 2009), storage calculation (Hollinger et al., 1994), and footprint analysis (Kljun et al., 2004) calculated using online Eddypro 4.2 software (Licor Inc., Lincoln, NE, USA). The Eddypro proceeding water vapor fluxes were flagged with different numbers, where 0 and 1 indicated high quality fluxes that were suitable for ET dynamics analysis, while 2 indicated poor quality fluxes that should be discarded. The Eddypro data for ET flagged as 0 and 1 were further filtered for precipitation, threshold deleting, and low turbulence. The ET fluxes selected during a precipitation event and the following half hour were deleted as the sensors may suffer from moisture contact. The selected ET lower than −0.32 mm h −1 or higher than 1.30 mm per hour were deleted. Low turbulence was determined by the u Ã threshold, which was calculated following Reichstein et al. (2005) using u Ã , CO 2 fluxes, and T a . The calculated u Ã thresholds were 0.25 and 0.21 m s −1 for 2016 and 2017, respectively. The data coverage after data quality control was 42.55% and 44.55% for 2016 and 2017, respectively. However, the data coverage showed obvious differences between daytime and nighttime. During the daytime, the data coverage was 66.66% and 66.28% for 2016 and 2017, respectively. However, nighttime only had 18.48% and 22.82% data coverage for 2016 and 2017, respectively. This was primarily affected by the low turbulence at night.
The ET data gaps were filled using the look-up table method (Falge et al., 2001a;Reichstein et al., 2005) based on T a , vapor pressure deficit (VPD), and R n .
CO 2 fluxes were measured simultaneously using the same EC system and were calculated as the covariance between vertical wind speed (v) fluctuations and CO 2 concentration fluctuations. Heat fluxes (H) were calculated as the covariance between vertical wind speed (v) fluctuations and air temperature fluctuations. The raw H fluxes suffered from quality control issues similar to those of ET.
G was also subjected to quality control by threshold values. Data gaps for G were linearly interpolated for gaps of less than 2 h but were filled using the mean diurnal variation (MDV) method for longer periods (Falge et al., 2001b).

Auxiliary data processing
The auxiliary data including conventional meteorological data and biotic data were also subjected to data quality controls and gap-fillings.
T a and RH above the canopy (at the height of 30 m) were used to analyze the drivers of ET dynamics. The qualities of T a and RH were firstly controlled through threshold values. The measured T a lower than −50 C or higher than 50 C were deleted, while the threshold values of RH were 0% and 100%. In addition, the measured T a and RH lower than 70% or higher than 130% of the mean corresponding values of other layers were deleted. The data gaps of meteorological variables (including T a and RH) were filled using the linear interpolation for gaps of less than 2 h. The remaining data gaps were firstly filled with the corresponding variables in neighboring layers. Then the remaining data gaps were filled with the mean diurnal variation (MDV) method (Falge et al., 2001b). VPD was then calculated as the difference between the actual and saturation vapor pressures based on the measured RH and T a .
Precipitation data were quality checked and any gaps were filled using manual observations from a meteorological station affiliated with the Chinese Bureau of Meteorology, which was located approximately 18 km away from the observation tower.
R g was quality controlled with spikes deletion. R g data lower than 0 or higher than 1,300 W m −2 were removed. The data gaps of R g were filled using the linear interpolation for gaps of less than 2 h, while the remaining data gaps were filled with the MDV method (Falge et al., 2001b).
R n was also quality checked with spikes deletion used to estimate the energy balance closure. First, short wave radiation data lower than 0 or higher than 1,300 W m −2 were removed. Second, long wave radiation data lower than 100 W m −2 or higher than 600 W m −2 were removed. Only data passing the quality control were used to calculate R n .
LAI was downloaded from the Moderate Resolution Imaging Spectroradiometer (MODIS) database with a spatial resolution of 500 m and a temporal variation of 4 days (https://modis.ornl.gov/subsetdata). In addition, the pixels from the surrounding 1 km area were also downloaded to avoid potential noises in a single pixel, which all had the same mixed forest. The downloaded LAI were linearly interpolated into a daily scale.

Energy balance closure assessment
We quantified the energy balance closure using two methods to validate the performance of the eddy covariance measurements. The first method derived the linear regression between the half-hour dependent flux variables (LE + H) and the available energy (R n − G) from the ordinary least squares (OLSs) relationship (Li et al., 2005). The regression coefficients, including the slope and the intercept, were used to evaluate the performance of energy balance closure. The ideal closure was reflected by a one slope and a zero intercept. The second method derived the cumulative sum of LE + H and available energy (R n − G) over measuring period and calculated the energy balance ratio (EBR) as follows: Only data directly measured and passing the quality control were used to calculate the energy balance closure. All available data were first potted to generate the energy balance closure during the measuring period. Then the whole year was divided into the active growing season (May-October) and the non-growing season (November to April), to generate the energy balance closure at different seasons.

ET dynamics and their drivers
Daily ET and environmental variables of each year were described to illustrate their seasonal dynamics. These were then used to reveal the drivers of the seasonal dynamics of ET. Only the data measured above the canopy (30 m) was used, although multilayer observations of T a , RH, and VPD were obtained. To avoid the confusion of gap-filling data in representing ET seasonal dynamics and their drivers, we only used daily data during days having more than 20 (including 20) directly measured and passing the quality control data. Setting 20 as the minimum data number sourced from the following four aspects. First, the missing of ET primarily resulting from low turbulence during nighttime made no day have all half-hour data directly measured. Second, data gaps were primarily caused by ET missing as nearly all R g were directly measured when the number of directly measured ET were over 20. Third, the small ET values during nighttime made the uncertainties in gap-filling contribute little to the daily accumulated values. Fourth, the occasionally appearing data gaps during daytime could be linear interpolated with high confidence.
The data from an entire year were divided into two groups: active growing season and non-growing season. The active growing season lasted from May to October and the non-growing season included the remaining months. We calculated the mean ET for each group over each half-hour and their corresponding variables, including T a , R g , and VPD, to represent the diurnal variations in different seasons. Only data directly measured and passing the quality control were used. Each available half-hour data point in each group was used to analyse the drivers of the diurnal variations of ET. Each year data were used seperately to validate the results from each year.
Though R n was the direct variable reflecting the available energy, the calculation of R n may be inaccurate as the bad sensors especially the long wave radiation sensors. Therefore, R g was used to replace R n to reveal the drivers of ET dynamics.

Uncertainty analysis
A simple Monte Carlo experiment was conducted to assess the uncertainty in the annual estimates of ET following Desai et al. (2005). A random number generator was used to remove 10-50% of the existing observed data. The generated data gaps were then randomly selected and filled with modeled ET (described in Flux processing), which was run 100 times for each year (2016 and 2017).

Statistical analysis
MATLAB 2014a (Mathworks Inc., Natick, MA, USA) was used to process the data including the flux quality control and gap filling. Linear regression and non-linear regression were used to analyse the effects of various factors on dynamics (including the diurnal and seasonal dynamics) of ET. Stepwise regression was used to analyze the multiple linear regression on ET seasonal variations. The minimum P-value for a variable to be added or removed from the model in the stepwise regression was 0.10. A path-analysis was conducted to evaluate the dependence of ET on various factors. The significant level was set to a = 0.05.

RESULTS
The energy balance closure The energy balance closure indicated by the regression slope between the available energy (R n − G) and energy fluxes (LE + H) (Fig. 2), suggested that our ecosystem experienced an energy imbalance. The slope between R n − G and LE + H during the measuring period was 0.54, with an intercept of 13.87 W m −2 ( Fig. 2A), which was obviously lower than 1. Additionally, different seasons had divergent regression slopes (Figs. 2B and 2C). During the active growing season, the regression slope could be 0.56 with an intercept of 21.68 W m −2 (Fig. 2B). Though the regression slope during active growing season was less than 1, it was still higher than that of the whole period (Fig. 3A) and the non-growing season, specifically (Fig. 2C). The regression slope during the non-growing season was only 0.39, with an intercept of 8.88 W m −2 (Fig. 2C).
The energy balance closure, as indicated by the energy balance ratio (EBR), also reflected the energy imbalance of the studied ecosystem. During the measuring period, EBR was 0.62 ( Fig. 2A), which was also clearly lower than 1. However, the EBR differed among seasons (Figs. 2B and 2C). During the active growing season, the EBR measured 0.66 (Fig. 2B), which was still lower than 1 but higher than that of the whole period ( Fig. 2A) and the non-growing season (Fig. 2C). The EBR during the non-growing season only measured 0.47 (Fig. 2C).

The seasonal dynamics of ET and environmental factors
ET and its related environmental factors all exhibited obvious seasonal dynamics (Fig. 3).
Environmental factors including T a , VPD, R g , precipitation, and LAI all exhibited obvious single peak seasonal dynamics (Figs. 3A-3F). T a ranged from −17 to 29 C, with the highest values appearing at July, which was consistent across both years measured (Fig. 3A). VPD had higher values during the active growing season and lower values during the non-growing season; the highest VPD appeared in May (Fig. 3C). R g ranged from 0.2 to 25 MJ m −2 d −1 and had similar dynamics to those of VPD (Fig. 3D). Precipitation also differed among seasons, with the most precipitation occurring from May to September (Fig. 3E). LAI ranged from 0.2 to 4.9 m 2 m −2 , with the highest values in July. RH did not show much variation across the months (Fig. 3B).
The seasonal dynamics of ET also exhibited a single-peak pattern, with the higher values appearing in the summer (Fig. 3F). However, the peak values of the seasonal dynamics of ET differed between the years measured (Fig. 3F). In 2016, ET achieved its peak (>6 mm d −1 ) in the middle of July and then decreased. In 2017, ET showed an increasing trend until the end of August, and achieved its highest value for 2017 (>7 mm d −1 ) at the end of August.

The diurnal dynamics of ET and environmental factors
ET and environmental factors all exhibited similar, obvious diurnal dynamics during the active growing season and non-growing season, but their values differed between seasons (Fig. 4).
During the day, T a and R g both showed single-peak patterns (Figs. 4A, 4B and 4G, 4H). These factors did not vary much from midnight to sunrise, but rapidly increased after sunrise, achieved their peak values around noon, and started to decrease at sunset. The factors were then stable until the end of the day. The time that T a achieved its peak values (approximately 13:00 local time) was later than that of R g . However, RH exhibited a concave shape (Figs. 4C and 4D), which indicated that RH began to decrease Figure 2 The energy balance closure calculated by the ordinary least squares (OLSs) relationship and energy balance ratio (EBR). The OLSs relationship was calculated between the available energy (R n -G) and the energy fluxes in the studied ecosystem during measuring period. The available energy was the difference between the net radiation (R n ) and soil heat flux (G). The energy fluxes were the sum of heat fluxes (H) and latent heat fluxes (LE). EBR was calculated as the ratio of total energy fluxes to available energy during the measuring period, was also calculated. Only data passing quality control were used. The energy balance closure was calculated during the whole measuring period (A), the active growing season (B), and non-growing season (C).
Full-size  DOI: 10.7717/peerj.13549/ fig-2 after sunrise. VPD did not vary much during the non-growing season but showed a single peak diurnal pattern in the active growing season (Figs. 4E and 4F). There were also similar diurnal dynamics for ET (Figs. 4I and 4J). ET showed single-peak diurnal patterns during both seasons, with the peak values occurring at 12:00 pm. However, the peak values differed between seasons. During the non-growing season, the ET peak value was no more than 0.05 mm h −1 , while the ET peak values during the active growing season could be 0.28 mm h −1 . These results were consisting across both years measured (Figs. 4I and 4J). Additionally, the highest diurnal peak value for ET differed between measured years, with the peak value during non-growing season in 2016 higher than that of 2017 (Fig. 4I), while the peak value during active growing season in 2016 was lower than that of 2017 (Fig. 4J).
In addition, the diurnal variations of ET and environmental factors showed some fluctuations (Fig. 4), which may source from the unequal number of available data in each half-hour.

Drivers of ET dynamics
Environmental factors including T a , R g , and VPD were found to significantly drive the single peak diurnal patterns of ET, but their effects differed among factors and seasons (Fig. 5). During the non-growing season, the increasing T a , R g , and VPD were found to significantly affect ET, while their effects differed among factors. T a exerted the strongest effects on ET diurnal variations during the non-growing season. The equations containing T a explained 32% and 26% of the variations for ET in 2016 and 2017, respectively (Figs. 5A and 5D). The effects of R g and VPD on the diurnal variations of ET during the non-growing season were similar but weaker than those of T a . Therefore, T a exerted a stronger effect on the diurnal variation of ET during the non-growing season. During the active growing season, the increasing ET was accompanied by an increase in T a , R g , and VPD, while their effects differed among factors. In 2016, T a exerted a stronger effect on ET diurnal variations during the active growing season of 2016 (Fig. 5G), while R g and VPD exerted weaker effects (Figs. 5H and 5I). The equation containing T a explained 24% of the diurnal ET variations. In 2017, R g had the strongest effect the diurnal variations of ET (Fig. 5J). The equation containing R g explained 42% of the diurnal ET variations. T a exerted a similar effect with R g , which were both stronger than VPD (Figs. 5J-5L). Therefore, the diurnal ET variations during the active growing season were governed by T a and R g , respectively.
The seasonal variation of ET was primarily shaped by T a and LAI, while R g and VPD did not contribute much (Fig. 6). ET increased exponentially as T a increased with an R 2 of 0.65 and 0.83 for 2016 and 2017, respectively (Fig. 6A). R g only significantly increased ET in 2017 with an R 2 of 0.15 (Fig. 6B).Though the increasing VPD significantly promoted ET, the equations containing VPD explained little of the ET seasonal variations (Fig. 6C). LAI made great contributions to the seasonal variations of ET, and as it increased, ET increased linearly with an R 2 of 0.69 for 2016 and 0.82 for 2017, respectively (Fig. 6D). The significant correlations between ET and the environmental variables ( Fig. 6) indicated that the seasonal variations of ET were shaped by multiple factors. Stepwise analysis showed that variables entering the equation describing ET seasonal variations differed between years. In 2016, only T a and LAI entered the equation describing ET seasonal variations, with an R 2 of 0.75 and an RMSE of 0.68 mm d −1 (Eq. (2)). In 2017, all variables including T a , R g , and LAI entered the equation describing ET seasonal variations, with an R 2 of 0.86 and an RMSE of 0.66 mm d −1 (Eq. (3)).
ET ¼ 0:038T a þ 0:604LAI þ 0:467; R 2 ¼ 0:75; RMSE ¼ 0:68 ET ¼ 0:030T a þ 0:031R g þ 0:754LAI À 0:092; R 2 ¼ 0:86; RMSE ¼ 0:66 Though many variables entered the equations describing the seasonal variations of ET (Eqs. (2) and (3)), each variable had a divergent role (Fig. 7). LAI exerted the strongest direct effect on the seasonal variation of ET. T a and R g also directly affected the seasonal variation of ET but showed a weaker effect than LAI. In addition, the seasonal variation of LAI was dominated by T a , and supported by R g and precipitation. The roles of environmental factors in the seasonal variation of ET were comparable in both years (Fig. 7). LAI was shown to exert a direct effect on the seasonal dynamics of ET in our ecosystem.

The annual amount of ET and its related variables
The mean annual ET during the measuring period was 528.26 mm year −1 , which were measured as 501.91 ± 5.30 mm year −1 and 554.60 ± 11.24 mm year −1 for 2016 and 2017, respectively (Table 1). The uncertainty of annual ET sourced from gap-filling was less than 2%. The active growing season values from May to October were the primary contributors to ET and accounted for over 90% of its annual value (Table 1). July had the highest ET, exceeding 120 mm month −1 in both years, which was 1 month prior to the measurement of the highest precipitation (Table 1). However, June contributed the second highest ET in 2016 while August had the second highest ET in 2017. The annual mean air temperature (MAT) were 0.57 and 1.17 C for 2016 and 2017, respectively, which resulted from a hot summer from May to October and a cold winter from November to April. The annual R g were not reported as many long data gaps during the measuring period, while VPD for 2016 and 2017 were only 0.33 and 0.37 kPa, respectively. The annual precipitation were 565 and 518.50 mm for 2016 and 2017, respectively. In addition, the annual mean LAI of 2016 and 2017 were 1.53 and 1.47 m 2 m −2 , respectively. Therefore, this forest was shown to experience a high evaporate rate, which is crucial to the global water cycle.

DISCUSSION
The quality of our measurements was key for the accurate analysis of ET dynamics and their drivers. These were evaluated using the available data portion and energy balance closure. The available data portion may be over 43.55% of the whole year as determined by the number of observations passing data quality control, which was slightly lower than that of some ChinaFLUX ecosystems (Fu et al., 2009;Yu et al., 2008;Zhu et al., 2015a). For example, the data coverage of three forests ranged from 43% to 54% (Zhu et al., 2015a), while that of three grasslands varied from 46% to 50% (Fu et al., 2009). Our ecosystem experienced a colder winter, which affected the eddy covariance instruments negatively, resulting in a lower data quality. The relative lower data coverage of our ecosystem was acceptable. From the energy balance closure, we found that the energy balance regression slope of our ecosystem from OLSs was 0.54 and the EBC was around 0.62, which was lower than most ChinaFLUX ecosystems (Li et al., 2005) and some European forests (Moderow et al., 2009;Sanchez, Caselles & Rubio, 2010) but fell into the range of energy balance closure of FLUXNET (Stoy et al., 2013;Wilson et al., 2002). The lower EBC in our ecosystem could be attributed to three aspects. First, our ecosystem experienced a cool temperature and a multiple plant composition, which made our ecosystem have a lower EBC as the EBC would increase with the increasing mean annual air temperature (Cui & Chui, 2019) and the decreasing heterogeneity (Foken, 2008;Stoy et al., 2013;Xin et al., 2018). Second, the soil heat flux (G) may be underestimated, especially during the non-growing season. This was also validated by an obviously higher energy balance closure during the active growing season (Figs. 2B and 2C). G was the product of soil bulk density, soil heat capacity, and soil temperature variations. Our ecosystem experienced a long snow covering period, whose heat capacity was higher than the soil (Li, Jia & Lu, 2015). However, we calculated the G during snow covering period with a soil heat capacity at ice free seasons, which made soil heat capacity used in calculating G during snow covering period underestimated. Therefore, the G during the non-growing season was underestimated, which further decreased the EBC of our ecosystem. Third, we only considered the heat storage in the soil (G) but ignored other heat storages. Heat can be stored in the soil, biomass, air (Moderow et al., 2009), and biogeochemical processes (Eshonkulov et al., 2019), with nearly equal contributions among these components (Lindroth, Mölder & Lagergren, 2010). Ignoring heat storage in some components (like air and biomass) may cause the energy balance closure to be underestimated (Eshonkulov et al., 2019;Franssen et al., 2010;Lindroth, Mölder & Lagergren, 2010). We can conservatively state that the eddy covariance measurement in our ecosystem performed well after fully considering the data coverage and energy balance closure. We analysed the dynamics of ET and their drivers based on EC measuring ET and environmental factors. Our results showed that the seasonal and diurnal variations of ET all exhibited a single peak pattern with the daily ET ranging from 0 to 7.75 mm d −1 and the hourly ET ranging from 0 to 0.28 mm h −1 . We also found that diurnal dynamics of ET during the non-growing season and the active growing season were driven by T a and R g , while seasonal ET dynamics were primarily affected by LAI. The single peak patterns of diurnal and seasonal dynamics for ET were commonly found in temperate ecosystems (Vourlitis et al., 2002;Wever, Flanagan & Carlson, 2002;Zheng et al., 2014;Zhou et al., 2010). However, our hourly and daily ET ranges differed from previous works, which may be related to a difference in climate and ecosystem types. For example, the daily ET of a Japanese temperate cypress forest ranged from 0 to 5.15 mm d −1 (Kosugi et al., 2007), while the daily ET of a Chinese warm temperate plantation ranged from 0 to 7.4 mm d −1 (Tong et al., 2017). American temperate managed forests had a daily ET ranging from 0 to 6 mm d −1 and hourly ET varying from 0 to 0.3 mm h −1 (Sun et al., 2008). The drivers also differed from other ecosystems (Tong et al., 2017;Yoshida et al., 2010), which may be related to the unique characteristics of our ecosystem. As a temperate mixed forest, our ecosystem experienced a unique climate with sufficient water supply but limited radiation. Therefore, energy supply, which could be represented by T a or R g , was the primary factor shaping the dynamics of ET in our ecosystem. T a showed a larger range than R g (Figs. 5A-5F), which made T a be the primary factor shaping ET diurnal dynamics when considering the diurnal variations of ET during the non-growing season. Given most data during active growing season of 2016 were missing (Fig. 3), the effects of environmental factors on ET diurnal variations were comparable to those during the non-growing season, which made T a be the dominating factor shaping ET diurnal variation. However, R g exhibited a larger range than T a (Figs. 5G-5L) during the active growing season of 2017, which indicated that R g governed the diurnal dynamics of ET. Though LAI was found to directly affect the seasonal dynamics of ET, LAI was the comprehensive representation of T a and R g (Fig. 7) and reflected the dominating role of energy supply in the seasonal dynamics of ET.
The unique climate resulted in an evaporation rate of 528.26 mm year −1 as ET, which was almost equal to its annual precipitation. The studied ecosystem had a lower ET than those with a warmer MAT like the Horqin grassland (Li et al., 2016). However, our ecosystem had a higher ET than other ecosystems in this region, including those experiencing a warmer MAT like the Changbaishan forest (Wu et al., 2015), Tomakomai forest (Hirano et al., 2003), and Laoshan forest (Wang et al., 2004) (Table 2), or ecosystems having a higher MAT and a lower MAP (including the Tongyu cropland, Tongyu grassland, Changling grassland) (Li et al., 2016;Liu & Feng, 2012), and ecosystems experiencing a similar MAT but a lower MAP like the Kherlenbayan-Ulaan grassland (Li et al., 2007;Liu et al., 2010). Our findings indicate that the studied ecosystem evaporated more water into the atmosphere, suggesting its importance in the global water cycle.

CONCLUSIONS
In this study, we analysed the dynamics of evapotranspiration (ET) in a temperate mixed forest using an eddy covariance approach. Our results showed that 43.55% eddy covariance measured data passed the data quality control checks with an energy balance ratio of 0.62. These results indicate the accuracy of the eddy covariance approach in our ecosystem. ET exhibited single-peak diurnal and seasonal patterns, with diurnal dynamics driven by air temperature (T a ) and global radiation (R g ) during non-growing season and active growing season, respectively, and seasonal dynamics affected by leaf area index (LAI), which all reflected the energy supply. The dynamics of ET resulted in a mean annual ET of 528.26 mm year −1 during 2016-2017. Therefore, the energy supply governed the dynamics of ET during all seasons and time scales in our temperate mixed forest, but variables representing the energy supply differed among seasons and time scales. The dynamics of ET show that this temperate ecosystem had an important role in global water cycles. The results of our work improve our understanding of ET dynamics in this region.