New gap-filling and partitioning technique for H₂O eddy fluxes measured over forests

The continuous measurement of H₂O fluxes using the eddy covariance (EC) technique is still challenging for forests because of large amounts of wet canopy evaporation (E_WC), which occur during and following rain events when the EC systems rarely work correctly. We propose a new gap-filling and partitioning technique for the H₂O fluxes: a model–statistics hybrid (MSH) method. It enables the recovery of the missing E_WC in the traditional gap-filling method and the partitioning of the evapotranspiration (ET) into transpiration and (wet canopy) evaporation. We tested and validated the new method using the data sets from two flux towers, which are located at forests in hilly and complex terrains. The MSH reasonably recovered the missing E_WC of 16–41 mm yr⁻¹ and separated it from the ET (14–23 % of the annual ET). Additionally, we illustrated certain advantages of the proposed technique which enable us to understand better how ET responds to environmental changes and how the water cycle is connected to the carbon cycle in a forest ecosystem.

1 Introduction

Forest ecosystems share three properties that are significant in their interactions with the atmosphere. They are extensive, dense, and tall, and thus produce sizable aerodynamic roughness and canopy storage for rainfall interception/evaporation (e.g., Shuttleworth, 1989). Since most flat terrain is used as agricultural or urban space, substantial areas of forest exist in mountainous terrains where the fundamental assumptions of eddy covariance (EC) measurement (flat and homogeneous site; e.g., Baldocchi et al., 1988) are invalidated. These facts hinder the use of the EC method from assessing the net ecosystem exchanges (NEE) of H₂O and CO₂ in forests.

Considering that EC measures compound “net” fluxes and that its gaps are unavoidable, we commonly take great care in flux gap-filling and partitioning. Basically, the gap-filling and partitioning are a kind of interpolation and extrapolation based on the fact that EC measurement has high temporal resolution and the bio-meteorological processes is a (repetitive) cycle (“redundancy” of data; Papale, 2012). Generally, they consist of the following procedure: (1) setting a target flux (e.g., CO₂ ∕ H₂O ∕ CH₄ fluxes, ecosystem respiration), (2) selecting drivers which control the target flux, (3) identifying relationships between the (appropriate) target flux (which can represent true NEE) and the drivers, and (4) interpolating and extrapolating the relationships during a certain period when the relationships are maintained (e.g., Papale, 2012; Reichstein et al., 2012). In this context, the gap-filling and partitioning (including nighttime CO₂ flux correction) are coterminous with each other. A related scientific issue is determining/selecting the number and type of drivers, and the method and the time window size used to identify the relationship. This depends on data availability, temporal scale of the process, and ecosystem state change. Those processes require extra care in measurement in complex mountainous terrain (e.g., van Gorsel et al., 2009; Kang et al., 2017).

Wet canopy evaporation (E_WC) is evaporation of the intercepted water by the vegetation canopy during and following rain events, which may consist of a significant portion of evapotranspiration (ET). Over forests, it is hard to measure the E_WC primarily due to the malfunction of an open-path EC system with rainfall. Although a closed-path system with an intake tube enables the E_WC measurement in the rain, the attenuation of the turbulent flow inside the tube acts as low-pass filtering, which results in a significant underestimation of the E_WC. Furthermore, the attenuation domain expands with increasing relative humidity (RH) from high to medium frequency (e.g., Ibrom et al., 2007; Fratini et al., 2012). The closed-path EC system with the heated tube may be the most appropriate for measuring ET in the rain (e.g., Goodrich et al., 2016).

The missing (or low quality) data can be gap-filled using general gap-filling methods such as marginal distribution sampling (MDS) and artificial neural network (e.g., Reichstein et al., 2005; Papale and Valentini, 2003). However, Kang et al. (2012) showed that, without proper consideration of the E_WC, such gap-filled ET data under the wet canopy conditions are underestimated because the data used in such gap-filling are mostly collected during dry or partially wet canopy conditions when the EC systems work properly. The authors proposed an improved gap-filing method that is coupled with a simple canopy (water) interception model.

The ET represents a combination of the E_WC, transpiration (T), and soil evaporation (E_S), which are controlled by different mechanisms and processes. Therefore, the partitioning of ET into the E_WC, T, and E_S is required to understand how ET is regulated by environmental changes and the how water cycle is connected to the carbon cycle in a forest ecosystem. For these reasons, there have been many studies that partition ET using other supplementary measurements or empirical/process models (e.g., Wilson et al., 2001; Yepez et al., 2003; Daikoku et al., 2008; Stoy et al., 2006; Hu et al., 2009; Kang et al., 2009b). Despite numerous previous studies on ET partitioning, most of them have focused on the partitioning of ET into the E_S (or direct evaporation, i.e., a sum of E_S and E_WC) and T. In the case of forest ecosystems with a dense canopy under a monsoon climate (e.g., East Asia, South Asia), E_WC can play a greater role than the E_S. In this context, it is necessary to pay attention to the method described by Kang et al. (2012), which not only allows for the proper estimation and gap-filling of the missing evaporation data under wet canopy conditions but also enables the partitioning of ET into the E_WC and T appropriately after certain modifications.

In this study, we propose a new gap-filling and partitioning technique for the H₂O fluxes measured over forests in complex mountainous terrain. First, we introduced a model–statistics hybrid (MSH) method, which can not only recover the missing E_WC in the general gap-filling method but also separate it from ET. Then, we tested and validated these new methods using the data sets from the two flux towers, which are located in forests with hilly and complex terrains. Additionally, we illustrated certain advantages of the new technique.

2 Materials and methods

2.1 Study sites

In the Korea National Arboretum, there are two eddy covariance flux towers: the Gwangneung deciduous forest located at the top of a hill (GDK; $\mathrm{37}{}^{\circ }{\mathrm{45}}^{\prime }{\mathrm{25}}^{\prime \prime }$ N, $\mathrm{127}{}^{\circ }{\mathrm{09}}^{\prime }{\mathrm{12}}^{\prime \prime }$ E) and the Gwangneung coniferous forest located at the bottom (GCK; $\mathrm{37}{}^{\circ }{\mathrm{44}}^{\prime }{\mathrm{54}}^{\prime \prime }$ N, $\mathrm{127}{}^{\circ }{\mathrm{09}}^{\prime }{\mathrm{45}}^{\prime \prime }$ E). Gwangneung has been protected to minimize human disturbance over the last 500 years. Both sites are located on complex, hilly catchment with a mean slope of 10–20^∘. The two towers are ∼ 1.2 km apart, and the mean slope between them is ∼ 6.2^∘ (Moon et al., 2005). The east/west slopes are gentle, whereas the north/south slopes are steep in the catchment. The mountain–valley circulation is the dominant wind regime in the sites (Hong et al., 2005; Yuan et al., 2007). Meteorological records from an automatic weather station ∼ 1.6 km northeast of the tower for 1997–2016 show that annual mean air temperature is 10.1 ± 0.6 ^∘C and the mean precipitation is 1472 ± 352 mm (National Climate Data Service System, http://sts.kma.go.kr/). At the GDK site, the vegetation is dominated by an old natural forest of Quercus sp. and Carpinus sp. (80–200 years old) with a mean canopy height of ∼ 18 m and a maximum leaf area index (LAI) of ∼ 6 m² m⁻² in June. Compared to the GDK site, the GCK site is in a lower area and is a flat, plantation forest with the dominant species of Abies holophylla (approximately 80 years old) with a mean canopy height of ∼ 23 m and a maximum LAI of ∼ 8 m² m⁻² in June. Further descriptions of the sites can be found in Kim et al. (2006) and Kang et al. (2017).

2.2 Measurements and data processing

The H₂O (and CO₂) fluxes have been measured since 2006 and 2007 at the GDK site and GCK site, respectively. At both sites, the EC system was used to measure the fluxes from a 40 m tower. The wind speed and temperature were measured with a three-dimensional sonic anemometer (model CSAT3, Campbell Scientific Inc., Logan, Utah, USA), while the H₂O (and CO₂) concentrations were measured with an open-path infrared gas analyzer (IRGA; model LI-7500, LI-COR, Inc., Lincoln, Nebraska, USA) at both sites. Half-hourly ECs and the associated statistics were calculated online from the 10 Hz raw data and stored in data loggers (model CR5000, Campbell Scientific Inc.). Other measurements such as net radiation, air temperature, humidity, and precipitation were sampled every second, averaged over 30 min, and logged in the data loggers (model CR3000 for the GDK site and CR1000 for the GCK site, Campbell Scientific Inc.). More information regarding the EC and meteorological measurements can be found in Kwon et al. (2009), and Kang et al. (2009a).

The multi-level profile systems were installed to measure the vertical profiles of the H₂O (and CO₂) concentrations at both sites and to estimate the storage flux using a closed-path IRGA (model LI-6262, LI-COR, Inc.). The measurement heights were 0.1, 1, 4, 8 (base of the crown), 12 (middle of the crown), 18 (the canopy top), 30, and 40 m for the GDK site and 0.1, 1, 4, 12 (base of the crown), 20 (middle of the crown), 23 (the canopy top), 30, and 40 m for the GCK site. More information regarding the multi-level profile system can be found in Hong et al. (2008) and Yoo et al. (2009).

To improve the data quality, the collected data were examined by the quality control procedure based on the KoFlux data processing protocol (Hong et al., 2009; Kang et al., 2014). This procedure includes a sector-wise planar fit rotation (PFR; Wilczak et al., 2001; Yuan et al., 2007, 2011), the WPL (Webb–Pearman–Leuning) correction (Webb et al., 1980), a storage term calculation (Papale et al., 2006), spike detection (Papale et al., 2006), gap-filling (MDS method; Reichstein et al., 2005), and nighttime CO₂ flux correction (van Gorsel et al., 2009; Kang et al., 2016). The details of the gap-filling and partitioning methods for the H₂O flux are described in the next chapters.

2.3 Gap-filling and partitioning methods for the H₂O flux

2.3.1 Marginal distribution sampling (MDS) method

The missing H₂O flux (i.e., evapotranspiration, ET) data were gap-filled using the MDS method (Reichstein et al., 2005; Hong et al., 2009). This method calculates a median ET under similar meteorological conditions within a time window of 14 days and replaces the missing values with the median. The intervals of the similar meteorological conditions were 50 W m⁻² for the downward shortwave radiation (R_sdn), 2.5 ^∘C for the air temperature (T_a), and 5.0 hPa for the vapor pressure deficit (VPD). If similar meteorological conditions were unavailable within the time window, its interval increased in increments of 7 days before and after the missing data point (i.e., 14-day window size) until it reached 56 days (i.e., before and after 7 days → 14 days → 21 days → 28 days). When the missing ET values could not be filled in a time window of less than 56 days, R_sdn was exclusively used following the same approach (i.e., calculating a median of ET under similar R_sdn conditions within a time window). This gap-filling method is used for not only the H₂O flux but also the sensible heat and daytime CO₂ fluxes.

2.3.2 Modeling of wet canopy evaporation

For estimating the wet canopy evaporation (E_WC), a simplified version of the Rutter sparse model (Valente et al., 1997) included in the VIC LSM (Variable Infiltration Capacity Land Surface Model; Liang et al., 1994) was used in the KoFlux data processing program. The E_WC is estimated as follows:

$$\begin{array}{}\text{(1)}& E_{\mathrm{WC}\mathrm{\_}\mathrm{Mod}}={\mathit{\sigma }}_{\mathrm{f}}{E}_{\mathrm{p}}{\left({\displaystyle \frac{W_{\mathrm{c}}}{S}}\right)}^n\left({\displaystyle \frac{r_{\mathrm{a}}}{r_{\mathrm{a}}+r_{\mathrm{0}}}}\right),\end{array}$$

where E_{WC_Mod} is the modeled E_WC, σ_f is the vegetation fraction (i.e., 1 minus the gap fraction); E_p is the potential evaporation ($E_{\mathrm{p}}=\frac{\mathit{\epsilon }A+\mathit{\rho }{c}_p\cdot \mathrm{VPD}\cdot g_{\mathrm{a}}/\mathit{\gamma }}{\mathit{\lambda }(\mathit{\epsilon }+\mathrm{1})}$, where ε is the dimensionless ratio of the slope of the saturation vapor pressure curve to the psychrometric constant γ, A is the available energy, ρ is the density of air, c_p is the specific heat of air, g_a is the aerodynamic conductance ($=\mathrm{1}/r_{\mathrm{a}})$, and λ is the latent heat of vaporization); r_a is the aerodynamic resistance to heat and water vapor transport; S is the canopy storage capacity; and r₀ is the architectural resistance. The term ($r_{\mathrm{a}}/(r_{\mathrm{a}}+r_{\mathrm{0}}))$ is added to consider the variation of the gradient of specific humidity between the leaves and the overlying air in the canopy layer. W_c is the intercepted canopy water, and the exponent n is an empirical coefficient.

W_c is estimated as

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial W_{\mathrm{c}}}{\partial t}}={\mathit{\sigma }}_{\mathrm{f}}P-D-E_{\mathrm{WC}\mathrm{\_}\mathrm{Mod}},\end{array}$$

where P is the input total rainfall and D is the drip. When W_c>0, the canopies are wet. When W_c>S, the drip starts (D>0).

Figure 1 Flowchart of the gap-filling and partitioning technique for evapotranspiration.

There are many inputs (i.e., E_p and P) and parameters (i.e., σ_f, S, n, r_a, and r₀) for estimating the E_WC and W_c. E_p, P, and $r_{\mathrm{a}}(=r_{\mathrm{am}}+r_b$ where r_am and r_b are the aerodynamic resistance of momentum transfer and the excess resistance; $r_{\mathrm{am}}=\stackrel{\mathrm{‾}}{U}/u_{\ast }^{\mathrm{2}}$ where $\stackrel{\mathrm{‾}}{U}$ is the wind speed, u^∗ is the friction velocity; and $r_b\approx \frac{\mathrm{4.63}}{u_{\ast }}$; Thom, 1972; Kim and Verma, 1990; Kang et al., 2009a) can be obtained/estimated from the flux tower measurement. The parameters can be divided into constant parameters (i.e., n and r₀) and seasonally varied parameters (i.e., σ_f and S). The default values (before optimization) of n and r₀ are two-thirds and 2 s m⁻¹, respectively. σ_f (=1 minus gap fraction) and S are functions of LAI (leaf area index): (1) the gap fraction is estimated by exp($-k\times$ LAI), where k varies from 0.3 to 1.5, depending on the species and canopy structure (Jones, 2013, k=0.75 and 0.485 for the GDK and GCK, respectively; (2) S is estimated by K_L × LAI, where K_L varies from 0.1 to 0.3 (default value of K_L=0.2; see Appendix A for more details). σ_f and LAI can be obtained from a plant canopy analyzer or digital photography (e.g., Macfarlane et al., 2007; Hwang et al., 2016). If actual measurement is not available, Moderate Resolution Imaging Spectroradiometer (MODIS) LAI can be used alternatively. In this study, σ_f (actually k) and LAI were estimated using a plant canopy analyzer (model LAI-2000, LI-COR, Inc.).

Table 1 Statistical parameters for the error assessment at the study sites. Mean bias error (MBE), mean absolute error (MAE), RMSE, and d indicate mean bias error, mean absolute error, root mean square error, and index of agreement, respectively. Slope and r² are from the linear regression analysis.

The generalization of the model can be augmented by providing the parameter optimization procedure using available flux data under wet canopy conditions. We argue that this is better than the validation using other data sets because the parameters may be site-specific (i.e., more validation does not fully guarantee the proposed model works properly everywhere). After optimizing the parameters (i.e., K_L, n, and r₀), the parameters changed slightly from the default values (see Appendix B for more details). Since the model results from before and after the parameter optimization were not statistically different in the error assessment, we still used the default values in a conservative way.

This method only considers the E_WC from the canopy by neglecting the E_WC from the trunk and stem. In addition, the interception of snow is not considered because the small amount of intercepted snowfall evaporates when the eddy covariance systems function improperly, and its melting and sublimation processes are much more complex than intercepted rainfall. To distinguish snowfall from total precipitation, the empirical discriminants in Matsuo et al. (1981) were used. This method uses air temperature and humidity near the ground surface to separate snow from rainfall because when it snows, air is not saturated and the near-ground air temperature is lower than that under rainy conditions. The result from this method should be scrutinized by comparing it with other precipitation data, which are measured at a weather station near the site.

2.3.3 Gap-filling and partitioning technique for evapotranspiration: MSH

The currently used MDS is expected to underestimate and overestimate ET under wet and dry canopy conditions, respectively, due to the gap-filling without the consideration of canopy wetness (because the evaporative fraction is proportional to canopy wetness). Therefore, the gap-filling technique for ET proposed by Kang et al. (2012) was used: (1) to identify the canopy wetness, the intercepted canopy water (W_c; see Eq. 1) was calculated using the simplified Rutter sparse model; (2) all the missing gaps were filled by the MDS using the data under dry canopy conditions only (i.e., when W_c=0), which corresponds to the ET under dry canopy conditions (ET_dry); (3) under wet canopy conditions (i.e., when W_c>0), the gap-filled data were replaced with the sum of the E_WC estimated by the simplified Rutter sparse model (i.e., E_{WC_Mod}) and the ET_dry multiplied by $\mathrm{1}-(W_{\mathrm{c}}/S{)}^n$ (i.e., the contribution from transpiration; see Eqs. 1 and 2).

Such a gap-filled ET was partitioned into the transpiration (T or ET from the dry canopy, which approaches the actual transpiration under a dense and closed canopy conditions) and E_WC as follows. When data were missing, the T was estimated as ($\mathrm{1}-(W_{\mathrm{c}}/S{)}^n)$ ET_dry, while the E_WC was estimated as E_{WC_Mod}. If the data were not missing (i.e., ET_Obs), the partitioning procedure was divided into two parts. If the signs of ET_Obs, ET_dry, and E_{WC_Mod} were the same, the T was estimated by multiplying ET_Obs and the ratios of ($\mathrm{1}-(W_{\mathrm{c}}/S{)}^n)$ ET_dry to the sum of ($\mathrm{1}-(W_{\mathrm{c}}/S{)}^n)$ ET_dry and E_{WC_Mod} (i.e., the estimated transpired fraction of ET), while the E_WC was estimated by multiplying ET_Obs and the ratios of E_{WC_Mod} to the sum of ($\mathrm{1}-(W_{\mathrm{c}}/S{)}^n)$ ET_dry and E_{WC_Mod} (i.e., the estimated evaporated fraction of ET). If the signs of ET_Obs, ET_dry, and E_{WC_Mod} were not the same, then T was estimated by ($\mathrm{1}-(W_{\mathrm{c}}/S{)}^n)$ ET_dry, while the E_WC was estimated by subtracting ($\mathrm{1}-(W_{\mathrm{c}}/S{)}^n)$ ET_dry (i.e., the estimated T) from ET_Obs. The procedure regarding the MSH is described in Fig. 1.

3 Results and discussion

3.1 Validation of the MSH

First, we evaluated the latent heat flux under (mostly) wet canopy conditions (λET_WC, i.e., λET when $W_{\mathrm{c}}/S>\mathrm{2}/\mathrm{3}$) from the MSH method (λET_{WC_MSH}) against the observed λET_WC (λET_{WC_Obs}) at both sites from 2008 to 2010 (Fig. 2). Since the λET_WC occurs in during and following rain events, open-path EC system can measure the λET_WC when the instrument dries out faster than the canopy. Most of the points are near the one-to-one line. The data scattered away from the one-to-one line are characterized by large aerodynamic conductance (e.g., > 100 mm s⁻¹) and/or large VPD (e.g., > 10 hPa).

Table 1 shows the statistical parameters for the error assessment (i.e., MBE, MAE, RMSE, d, slope, and r²; see Appendix D for more details about the error assessment). The slopes from the linear regression analysis are 0.97 ± 0.15 and 0.89 ± 0.07 with 0.69 ± 0.06 and 0.81 ± 0.02 of r² for the GDK and GCK sites, respectively. The d values for the sites were close to 1 (0.91 ± 0.01 for the GDK and 0.95 ± 0.01 for the GCK). Compared to previous research (i.e., Kang et al., 2012), the results from the MSH were closer to the observation due to the consideration of ET from the dry canopy. One of the leading causes of the error in λET_{WC_MSH} was identified as the discrepancy between the time when the rain occurred and the time the tipping bucket was tipped. The results from the further evaluation of the MSH using a closed-path EC system were similar to those using the open-path EC system (see the Appendix C). To validate only E_WC, cross-validation using the other models (e.g., Gash sparse analytical model; Gash et al., 1995) can be attempted (e.g., Kang et al., 2012). Overall, the results from the linear regression analysis of λET_{WC_MSH} and λET_{WC_Obs} show that MSH can provide λET_WC reasonably well for the sites.

Figure 2 Comparison of the latent heat flux under (mostly) wet canopy conditions (i.e., $W_{\mathrm{c}}/S>\mathrm{2}/\mathrm{3}$, where W_c is the intercepted canopy water and S is the canopy storage capacity) at the GDK (a) and GCK (b) sites: λET_{WC_Obs} indicates the observed latent heat flux under wet canopy conditions (λET_WC), while λET_{WC_MSH} indicates the estimated λET_WC using the MSH. The dotted line represents the 1:1 line.

3.2 Comparison between the MDS and the MSH

To evaluate the superiority of the MSH, we filled in the missing λET_WC data by using the MDS (λET_{WC_MDS}) and the MSH (λET_{WC_MSH}). The underestimation of the λET_{WC_MDS} was shown by the comparison with the sum of energy flux components except for latent heat flux (= net radiation + sensible heat flux + storage flux) in our previous study (Kang et al., 2012). The λE_{WC_mod} displayed the mirrored patterns of the sum of the other energy budget components, while the λET_{WC_MDS} was very small (mainly due to the low radiation during the rainy days). Thus, we expected that the MDS underestimates the ET since it cannot explicitly consider the key processes of wet canopy evaporation (i.e., the effects of aerodynamic conductance (g_a) change and sensible heat advection; see Kang et al., 2012 for a more detailed explanation). Actually, the average annual MBEs from 2008 to 2010 were −18 ± 6 W m⁻² for the GDK site and −15 ± 5 W m⁻² for the GCK site, respectively. It also should be noted that some λET_{WC_MSH} varied while λET_{WC_MDS} was nearly constant, because (1) the λET_{WC_Obs} rarely existed close to the missing data and (2) the MDS did not consider the effect of g_a (not shown here).

Figure 3 Seasonal variation of monthly integrated evapotranspiration (ET) with the gap-filled by the MDS method (ET_MDS); the ET gap-filled by the MSH method (ET_MSH), transpiration and wet canopy evaporation partitioned by the MSH method (T_MSH and E_{WC_MSH}), for the GDK (a) and GCK (b) sites. E_{S_Obs} indicates soil evaporation measured by the supplementary eddy covariance systems at the floors (adapted from Kang et al., 2009b).

Figure 4 Wavelet coherence spectrum of evapotranspiration (ET) with rainfall (P) for the GDK site. The thick solid contour is the 5 % significance level against red noise as calculated from a Monte Carlo simulation. Arrows are the relative phase angle with in-phase angle (positive correlation) pointing right, anti-phase angle (negative correlation) pointing left, and P leading ET by 90^∘ pointing down. The shaded area indicates the cone of influence where the edge effects might distort the results.

Figure 3 shows the monthly ETs gap-filled by the MDS and MSH methods for the GDK and GCK sites. First, the annual ETs from the MSH method were 16–41 mm yr⁻¹ larger than those from the MDS method, while the random uncertainties in gap-filled annual ETs were approximately 5 mm yr⁻¹ for both sites (quantified according to Finkelstein and Sims, 2001, and Richardson and Hollinger, 2007). Significant differences were identified in June, July, August, and September during intense rainfall. The biggest difference is shown in 2010 with more frequent and larger rainfall (for the GDK, the number of rainy days is 86, 82, and 103 days and the total amount of rainfall is 1407, 1323, and 1652 mm in 2008, 2009, and 2010, respectively; such characteristics are similar to those for the GCK). In addition to taking the missing E_WC from the MDS into account, the other advantage of the MSH method is that the observed ET using the eddy covariance system can be partitioned into transpiration (T) and E_WC without any additional measurement. However, it can be applied to a dense canopy only, where soil evaporation is negligible. Otherwise (e.g., before leaf unfolding and after leaf fall), the T includes the error of the soil evaporation (E_S). Thus, there is more separating the E_WC than partitioning the ET. The annual E_WC ranged from 53 to 82 mm for the GDK and 78 to 112 mm for the GCK, which occupies 14–23 and 14–19 % of the annual ET, respectively.

For quantifying the E_S, the supplementary eddy covariance (EC) systems were operated at the floors of the GDK and GCK sites (Kang et al., 2009b). The annual understory ET (∼ E_S) from 1 June 2008 to 31 May 2009 was 59 mm for the GDK and 43 mm for the GCK, which comprised 16 and 8 % of the annual ET, respectively. The decoupling factor (Ω, McNaughton and Jarvis, 1983) at the forest floor was ∼ 0.15 for both sites, which indicates that the E_S was controlled primarily by the VPD and surface conductance (g_s) rather than R_sdn. This factor also suggests that separating E_S from ET using the exponential radiation extinction model to estimate the R_sdn at the forest floor and the relationship between the estimated R_sdn and the ET when the canopy is inactive (Stoy et al., 2006) can be problematic for the sites. Considering that the accurate estimation of g_s is challenging, a supplementary measurement (e.g., low-level EC, lysimeter, sap-flow and isotope measurements) is a better approach for estimating E_S. Using the E_S measured by the low-level EC, the annual T can be estimated at 265 mm (70 % of ET) for the GDK and 448 mm (78 % of ET) for the GCK, while the E_WC is estimated as 55 mm (15 % of ET) and 82 mm (14 % of ET), respectively (Fig. 3).

Table 2 The annual CO₂ and H₂O budget (NEE, net ecosystem exchange; GPP, gross primary production; RE, ecosystem respiration; ET, evapotranspiration; E_WC, wet canopy evaporation) and water use efficiency at the ecosystem level (WUE_Eco) and the canopy level (WUE_Canopy) for the study sites.

In the following chapters, we illustrate the advantages of the proposed technique. The benefits are caused by the gap-filling and partitioning of the H₂O flux because the MSH method can take the E_WC into account properly and separate them from ET, which has not yet been previously possible. We hope the following chapters draw attention to the ET partitioning.

3.3 Application: wavelet coherence analysis between ET and the rainfall

To evaluate the effect of new gap-filling, we conducted the wavelet coherence analysis between ET and rainfall for the GDK site (Fig. 4; see Hong et al., 2011, and Grinsted et al., 2004, for more details regarding the wavelet coherence analysis). From a 1- to 3-month period during the monsoon season (i.e., the intense rainy period), high correlation (i.e., red color area) was observed between ET and rainfall in 2006, 2007, 2008, and 2009. In 2007, the rainfall amount was 200 mm lower than the average level during the study period. However, the rainfall duration was the longest, and the intensity was the lowest. In 2006, 2008, and 2009, the arrow on the high correlation area pointed left. This means a negative correlation between the two variables, reflecting that the decrease in T was caused by the diminishment of R_sdn during the intense rainy period. In contrast, the arrow pointed right in 2007, indicating a positive correlation. The magnitude of enhanced E_WC was greater than that of decreased T at that time and frequency (i.e., August and September, 1- to 3-month period) in 2007. Such a positive correlation between ET and rainfall with the 1- to 3-month cycle in 2007 was not reported in the study of Hong et al. (2011), which showed a negative correlation in 2006, 2007, and 2008 at that time and frequency. This can be attributed to the improvement in ET data made by the new gap-filling method (i.e., recovering the missing E_WC in the general gap-filling method). During the monsoon season, the E_WC compensates for (a portion of) the decreased T, and are occasionally be balanced (e.g., in 2010).

There is some circumstantial evidence which suggesting that the proposed method is more appropriate for taking E_WC into account than the conventional method: (1) the ratio of the runoff and the precipitation (adapted from Choi, 2011) in 2007 was the lowest (0.60 in 2007, 0.69 ± 0.06 in the other years; i.e., the ratio of ET to precipitation was highest in 2007), while the R_sdn (main controlling factor of T) was lowest (4.52 GJ m⁻² in 2007, 4.77 ± 0.08 GJ m⁻² in the other years) due to the longest rainfall duration, (2) the interannual variabilities of the estimated catchment-scale annual ET (i.e., precipitation–runoff) and ET from the MDS method occurred in opposite directions (similarly to T from the MSH method).

3.4 Application: water use efficiency at the ecosystem level and the canopy level

Water use efficiency (WUE) can be defined in various forms such as A_n∕g_st (intrinsic WUE; A_n: net assimilation; g_st: stomatal conductance; A_n= NPP), A_n∕T (instantaneous WUE), $A_{\mathrm{n}}(\mathrm{1}-{\mathrm{\Phi }}_{\mathrm{c}})/[T$(1+Φ_w)] (integrated WUE; Φ_c: fraction of assimilated carbon lost in respiration; Φ_w: fraction of total water loss from non-photosynthetic parts of the plant or through open stomata at night), the GPP ∕ T (canopy-level WUE), NPP ∕ ET (stand-level WUE), and GPP ∕ ET (ecosystem-level WUE; Seibt et al., 2008; Ito and Inatomi, 2012; Ponton et al., 2006), because the spatiotemporal scale and measurement method are research-specific. Based on the original definition of WUE (i.e., the ratio of CO₂ flux to H₂O flux), we redefined the annual ecosystem-level WUE (WUE_Eco) and the annual canopy-level WUE (WUE_Canopy) as ΣNEP ∕ ΣET and ΣNPP ∕ ΣT, respectively. For estimating ΣNPP and ΣT simply, we used 0.45 as the ratio of the NPP to GPP for both sites (Waring et al., 1998), and 0.156 and 0.075 as the ratios of E_S to ET for the GDK and GCK, respectively (Kang et al., 2009b). From 2006 to 2010, WUE_Eco (WUE_Canopy) ranged from −0.16 (2.17) to 0.32 (2.59) g C (kg H₂O)⁻¹ for the GDK site and from 0.20 (1.93) to 0.38 (2.16) g C (kg H₂O)⁻¹ for the GCK site (Table 2). Considering the increasing trend of NEE and GPP for the GCK site, we identified that the interannual variabilities of WUE_Eco and WUE_Canopy occurred in opposite directions for both sites. It was primarily caused by E_WC being enhanced in 2007 and 2010 due to the weakest rainfall intensity and the largest rainfall amount, respectively. Overall, such partitioning of the total ET into E_WC, T, and E_S enables us to understand better how ET responds to environmental changes and the how water cycle is connected to the carbon cycle in a forest ecosystem.

4 Conclusions

The new technique proposed in this study for gap-filling and partitioning of H₂O eddy fluxes has a special feature: two existing methods were merged into a new method. The marginal distribution sampling (MDS) method and the simplified Rutter spars model have been merged into the model–statistics hybrid (MSH) method. Such a strategy increases the strength of the original methods while making up for their weaknesses. Thus, further improvement is necessary and forthcoming. The MSH method can be applied to tropical forests because tropical forests also share three properties of temperate forests (i.e., extensive, dense, and tall). However, applying the methods to grasslands may need further validation.

Appendix A: Parameterizations of canopy storage

The rainfall interception is sensitive to the change of canopy storage capacity (S) and vegetation fraction (σ_f, i.e., 1 minus the gap fraction; e.g., Shi et al., 2010). One of the characteristics common among the temperate forests considered in this study is a dense canopy. It means that the vegetation fractions in the forests are close to 1 (except before leaf unfolding and after leaf fall periods). Moreover, the gap fraction can be measured using a plant canopy analyzer with relative ease. Therefore, the error of the result from the model is mostly derived from the parameterization of S (see Appendix B for more detailed information). The canopy storage capacity is affected by not only leaf area but also the other factors such as leaf shape, leaf angle, leaf/shoot clumping, and hydrophobicity (water repellency) of a leaf (e.g., Crockford and Richardson, 2000). Additionally, the relationships between S and these characteristics are changeable according to meteorological conditions (e.g., the wind, rainfall intensity), which make the parameterization of S difficult (e.g., Dunkerley, 2009). Therefore, we used the simple parameterization of S of VIC LSM (i.e., $S=K_L\times$ leaf (or plant) area index, where K_L=0.2; Liang et al., 1994) and evaluated whether the parameterization was reasonable or not. The relationships between the leaf area index (LAI) and S in previous studies are presented in Fig. A1, indicating that the parameterization in VIC LSM (i.e., K_L=0.2) is reasonable, and the K_L ranges from 0.1 to 0.3. Further studies on the parameterization of S using leaf structure (e.g., leaf shape, leaf angle, leaf/shoot clumping) would be worth conducting for more accurate estimation of wet canopy evaporation.

Figure A1 Relationship between canopy storage capacity and plant/leaf area index (the data obtained from Table A1).

Table A1 Review of the canopy storage capacity for wet canopy evaporation modeling in previous studies.

Appendix B: Sensitivity test and parameter optimization of the wet canopy evaporation model

In the simplified Rutter sparse model (Liang et al., 1994), there are many parameters (e.g., σ_f, S, n, r_a, and r₀) for estimating the wet canopy evaporation (E_WC) and the intercepted canopy water (W_c). Since the model results may be sensitive to the parameters, and the parameters may be site-specific, the parameter optimization using available flux data under wet canopy conditions should be accompanied by the model for the generalization of the model. Considering that the gap-filling and partitioning are a kind of interpolation and extrapolation (i.e., identifying relationships between a target flux and its drivers, and interpolating and extrapolating the relationships), it is an appropriate strategy for the gap-filling and partitioning of evapotranspiration using the model.

First, we conducted a sensitivity test of the model to the parameters (i.e., k, K_L, n, and r₀) using the data set in 2008 (change in E_WC (%) $=\frac{E_{\mathrm{WC}\mathrm{\_}\mathrm{perturb}}-E_{\mathrm{WC}\mathrm{\_}\mathrm{default}}}{E_{\mathrm{WC}\mathrm{\_}\mathrm{default}}}\times \mathrm{100}$; E_{WC_default}: annually integrated E_WC simulated with default parameters; E_{WC_perturb}: annually integrated E_WC simulated after a change in each parameter. Only one parameter is changed one at a time, while other parameters are held constant; e.g., Shi et al., 2010). Before testing the sensitivity, we set the lower/upper boundaries (and default values) based on literature reviews: k=0.3 ∼ 1.5 (Jones, 2013; the default values of k are 0.75 and 0.485 for the GDK and GCK, respectively; these values were obtained from the actual measurement using a plant canopy analyzer, LAI-2000 from LI-COR, Inc.); K_L=0.1 ∼ 0.3 (see Appendix B; default value of K_L is 0.2, Dickinson, 1984); n=0.5 ∼ 1 (Chen and Dudhia, 2001; Liang et al., 1994; Valente et al., 1997; the default value of n is two-thirds: Deardorff, 1978); r₀= 0 (for short vegetation) and ∼ 2 s m⁻¹ (for tall vegetation; Perrier, 1975; Rana et al., 1994; the default value of r₀ is 2: Perrier, 1975). K_L is the most influential parameter (Fig. B1), implying that we should take great care to minimize its parameter estimation error. Using a small number of the observed latent heat flux data under wet canopy conditions (when $W_{\mathrm{c}}/S>\mathrm{2}/\mathrm{3}$) from 2008 to 2010, we optimized the parameters except k (because we obtained the k from the actual measurement) towards minimizing the root mean square error of the method (using the bound constrained optimization code in MATLAB^®, “fminsearchbnd.” http://kr.mathworks.com/matlabcentral/fileexchange/8277-fminsearchbnd–fminsearchcon). We randomly divided the available data set into the data sets for parameter optimization and validation (i.e., validation after optimization). The ratio of the optimization–validation data sets was arbitrarily set to 7:3. Table B1 shows the model parameters and the statistical parameters for the error assessment before and after the parameter optimization. After the optimization, the parameters were slightly different from the default values. However, we still used the default values conservatively since the model results from before and after the optimization were not statistically different in the error assessment.

Figure B1 Sensitivity test of the wet canopy evaporation model to the parameters (i.e., k, K_L, n, and r₀).

Table B1 Statistical parameters before and after the parameter optimization for the error assessment at the study sites. MBE, MAE, RMSE, and d indicate mean bias error, mean absolute error, root mean square error, and index of agreement, respectively. Slope and r² are from the linear regression analysis. The default values of K_L, n, and r₀ were 0.2 (0.2), 2∕3 (2∕3), and 2 (2) for the GDK (GCK), respectively. After the optimization, those values changed to 0.1966 (0.2314), 0.7279 (0.6930), and 2 (2) for the GDK (GCK), respectively.

Appendix C: Further evaluation of the MSH using a closed-path EC system

For further evaluation of the MSH method, we conducted a commercialized closed-path infrared gas analyzer with the heated tube (model EC155, Campbell Scientific Inc., Logan, Utah, USA) from 18 August 2015 to 25 October 2015 for the GDK site. The open- and closed-path gas analyzers shared the one sonic anemometer. We applied the MSH to the latent heat fluxes (λET) from the open-path eddy covariance (EC) system for gap-filling the λET under wet canopy conditions, and the frequency response correction (Fratini et al., 2012) to the λET from the closed-path EC system for correcting the tube attenuation effect especially under high relative humidity (RH) conditions. Then, we compared the λET (mostly) wet canopy conditions (λET_WC, i.e., λET when $W_{\mathrm{c}}/S>\mathrm{2}/\mathrm{3}$) from the MSH (λET_{WC_MSH_OP}) against the observed λET_WC from the closed-path EC system (λET_{WC_Obs_CP}) similar to Sect. 3.1 (Fig. C1).

Before the comparison, it should be noted that the data retrieval rate under wet canopy conditions of the closed-path EC system was 50 % higher than that of the open-path EC system; however, there were still data missing for a considerable period despite the use of the closed-path gas analyzer. The lack of data was mainly caused by the malfunction of the sonic anemometer and the unsatisfactory conditions for EC measurement (e.g., nonstationary and unfavorable turbulent conditions developed). The results of the error assessment (i.e., 10 W m⁻² of MBE, 19 W m⁻² of MAE, 29 W m⁻² of RMSE, 0.90 of d, 1.04 of slope, 0.68 of r²) were within the ranges of those from the comparison between the λET_{WC_MSH_OP} and the λET_{WC_Obs} in Sect. 3.1 (i.e., the case of the open-path EC system). Overall, such results imply the robustness of the MSH method as well as the necessity of appropriate λET_WC gap-filling method (e.g., MSH) when measuring using a closed-path EC system, as in the case of an open-path EC system.

Figure C1 Comparison of the latent heat flux under (mostly) wet canopy conditions (i.e., $W_{\mathrm{c}}/S>\mathrm{2}/\mathrm{3}$, where W_c is the intercepted canopy water and S is the canopy storage capacity) at the GDK site: λET_{WC_Obs_CP} indicates the observed latent heat flux under wet canopy conditions (λET_WC) from a closed-path eddy covariance (EC) system, while λET_{WC_MSH_OP} indicates the estimated λET_WC from an open-path EC system using the MSH method. The dotted line represents the 1:1 line.

Appendix D: Error assessment

In order to evaluate the latent heat flux under wet canopy conditions obtained from the MSH method, we compared it against the observed data using four statistical measures, following Willmott and Matsuura (2005). Mean bias error is the average of the residuals. Mean absolute error is the average of the absolute values of the residuals. A large deviation from zero implies that the estimation generally overestimates or underestimates compared to the observed values. We also considered root mean squared error (RMSE) which is often reported with MAE because RMSE is more sensitive to large errors than MAE.

$$\begin{array}{}\text{(D1)}& {\displaystyle }& {\displaystyle }\mathrm{MBE}=\sum {\displaystyle \frac{Y_{\mathrm{est}}-Y_{\mathrm{obs}}}{n}}\text{(D2)}& {\displaystyle }& {\displaystyle }\mathrm{MAE}=\sum {\displaystyle \frac{\left|Y_{\mathrm{est}}-Y_{\mathrm{obs}}\right|}{n}}\text{(D3)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{\sum {\displaystyle \frac{(Y_{\mathrm{est}}-Y_{\mathrm{obs}}{)}^{\mathrm{2}}}{n}}}\end{array}$$

MBE, MAE, and RMSE give estimates of the average error, but none of them provide information about the relative size of the average difference. Thus, we further considered an additional index of agreements (d), following Willmott (1982):

$$\begin{array}{}\text{(D4)}& {\displaystyle }d=\mathrm{1}-\left[{\displaystyle \frac{\sum (Y_{\mathrm{est}}-Y_{\mathrm{obs}}{)}^{\mathrm{2}}}{\sum (\left|Y_{\mathrm{est}}^{\prime }\right|+\left|Y_{\mathrm{obs}}^{\prime }\right|{)}^{\mathrm{2}}}}\right],\end{array}$$

where $Y_{\mathrm{est}}^{\prime }=Y_{\mathrm{est}}-\stackrel{\mathrm{‾}}{Y_{\mathrm{obs}}}$ and $Y_{\mathrm{obs}}^{\prime }=Y_{\mathrm{obs}}-\stackrel{\mathrm{‾}}{Y_{\mathrm{obs}}}$ (where overbar is an averaging operator). It ranges from 0 to 1, where 0 is for complete disagreement and 1 for complete agreement between the observation and the estimates. It is both a relative and bounded measure that can be widely applied in order to make a cross-comparison between models.