New techniques for gap-filling and partitioning of H 2 O and CO 2 eddy fluxes measured over forests in complex mountainous terrain

The continuous measurement of H2O and CO2 fluxes using the eddy covariance (EC) technique is still challenging for forests in complex terrain because of large amounts of wet canopy evaporation (EWC), which occur during and following rain events when the EC systems rarely work correctly, and the horizontal advection of CO2 generated at night. We propose 15 new techniques for gap-filling and partitioning of the H2O and CO2 fluxes: (1) a model-stats hybrid method (MSH) and (2) a modified moving point test method (MPTm). The former enables the recovery of the missing EWC in the traditional gap-filling method and the partitioning of the evapotranspiration (ET) into transpiration and (wet canopy) evaporation. The latter determines the friction velocity (u*) threshold based on an iterative approach using moving windows for both time and u*, thereby allowing not only the nighttime CO2 flux correction and partitioning but also the assessment of the significance of the 20 CO2 drainage. We tested and validated these new methods using the datasets from two flux towers, which are located at forests in hilly and complex terrains. The MSH reasonably recovered the missing EWC of 16 ~ 41 mm year-1 and separated it from the ET (14 ~ 23% of the annual ET). The MPTm produced consistent carbon budgets using those from the previous research and diameter increment, while it has improved applicability. Additionally, we illustrated certain advantages of the proposed techniques, which enables us to understand better how ET responses to environmental changes and how the water cycle is 25 connected to the carbon cycle in a forest ecosystem.


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., 30 Shuttleworth, 1989).Since most of the flat terrains are used as agricultural lands and towns (or cities), substantial areas of Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.
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 violated.These facts hinder the use of the EC method from assessing the net ecosystem exchanges (NEE) of H2O and CO2 in forests.

Wet canopy evaporation: gap-filling and partitioning
Wet canopy evaporation (EWC) is an evaporation of the intercepted water by the vegetation canopy during and following rain 5 events, which may consist of a significant portion of evapotranspiration (ET).Over forests, it is hard to measure the EWC primarily due to the malfunction of an open-path EC system with rainfall.Although a closed-path system with an intake tube enables the EWC measurement in the rain, the attenuation of the turbulent flow inside the tube acts as a low pass filtering, which results in a significant underestimation of the EWC.Furthermore, the attenuation domain expands with an increasing relative humidity (RH) from high frequency to medium frequency (e.g., Ibrom et al., 2007;Fratini et al., 2012).10 The missing (or low quality) data can be gap-filled using general gap-filling methods such as the marginal distribution sampling and artificial neural network (e.g., Reichstein et al., 2005;Papale and Valentini 2003).However, Kang et al. (2012) showed that, without a proper consideration of the EWC, such gap-filled ET data under the wet canopy conditions are underestimated because the data used in such gap-filling are mostly collected dry or partially wet canopy condition when the EC systems work properly.The authors proposed an improved gap-filing method that is coupled with a simple canopy (water) interception model.15 The ET represents a combination of the EWC, transpiration (T), and soil evaporation (ES), which are controlled by different mechanisms and processes.Therefore, the partitioning of ET into the EWC, T, and ES 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 previous 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., 2009a).20 Despite the many previous studies on ET partitioning, most of them have focused on the partitioning of ET into the ES (or direct evaporation, i.e., a sum of ES and EWC) and T. In the case of forest ecosystems with a dense canopy under a monsoon climate (e.g., East Asia, South Asia), the EWC can play a greater role than the ES.In this context, it is necessary to pay attention to the method described by Kang et al. (2012), which not only allows the proper estimation and gap-filling of the missing evaporation data under wet canopy conditions but also enables the partitioning of ET into the EWC and T appropriately after 25 certain modifications.

Nighttime CO2 flux: correction for complex topography
In EC measurements over complex mountainous terrain, the nighttime CO2 flux correction is one of the most important and challenging tasks.There are two types of widely used methods: (1) the friction velocity (u * ) filtering method, and (2) the advection-based filtering method.The most commonly used method is the u * filtering method that optimizes the parameters 30 of the ecosystem respiration (RE) function using the observed nighttime CO2 flux when u * is greater than a threshold, (i.e., there is no dependency of the nighttime CO2 flux on u * ) (e.g., Falge et al., 2001;Gu et al., 2005).The filtered data are replaced Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.
with the estimated data using ecosystem temperatures and the RE function with the optimized parameters.The u * filtering method cannot be applied at sites where the u * threshold cannot be determined and/or the drainage flow is generated at night.
Accordingly, an advection-based method was developed for hilly terrain sites that are affected by drainage flow (van Gorsel et al., 2007(van Gorsel et al., , 2008(van Gorsel et al., and 2009)).It uses the observed CO2 flux data near sunset when the nighttime advection has not yet affected to EC measurement.The nighttime correction is also important for the partitioning of NEE into gross primary productivity 5 (GPP) and RE because the nighttime RE-temperature relationship is used to extrapolate the daytime RE (e.g., Reichstein et al., 2005).
The Gwangneung deciduous and coniferous forests sites in Korea, which are a part of Korea Flux Monitoring Network (KoFlux), are typical sites situated in hilly and complex terrains where these methods are difficult to apply appropriately because the CO2 drainage generated earlier than the time assumed by the method.Kang et al. (2017) developed the site-specific 10 quality control filter to exclude the data strongly affected by CO2 advection.The filter identifies the observations that occur when a strong information flow toward the bottom of the slope exists in the dynamical process network of the observed multilevel CO2 concentrations.This site-specific filter, which is qualitatively similar to the application of the traditional correction methods in a hybrid way, substantially reduced the disagreement among the three different conventional methods of nighttime CO2 corrections.However, this method has low applicability because the measurements of CO2 concentration profiles are 15 required at two or more locations along the drainage, and a long time series data is necessary to produce robust results.To overcome these shortcomings, we propose a hybrid of the u * filtering and advection-based methods, which can be applied unconditionally, 'everywhere, all of the time.'

Overview of research
In this study, we propose new techniques for gap-filling and partitioning of the H2O and CO2 fluxes measured over forests in 20 complex mountainous terrain.First, we introduce a model-stats hybrid method, which can not only recover the missing EWC in the general gap-filling method but also separate it from ET.Then, an automated statistical method is introduced to determine the u * threshold based on an iterative approach using a moving window of both time and u * (i.e., the modified moving point test method).This method enables the determination of 'the significance of CO2 drainage' and the u * threshold for the nighttime CO2 flux correction and partitioning.We tested and validated these new methods using the datasets from the two flux towers, 25 which are located in forests with hilly and complex terrains.Additionally, we illustrated certain advantages of the new techniques.

4
37° 44' 54" N, 127° 09' 45" 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 are ~ 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 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 5 mean air temperature is 10.1±0.6°C and the mean precipitation is 1,472±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 2 m -2 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 2 m -2 in June. 10 Further descriptions of the sites can be found in Kim et al. (2006) and Kang et al. (2017).

Measurements and data processing
The H2O and CO2 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 H2O and CO2 15 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 dataloggers (Model CR5000,Campbell Scientific Inc.).Other measurements such as net radiation, air temperature, humidity, and precipitation were sampled every second, averaged over 30 minutes, and logged in the dataloggers (Model CR3000 for the GDK site and CR1000 for the GCK site, Campbell Scientific Inc.).More information regarding the EC and 20 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 CO2 and H2O 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 25 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;Yuan et al., 2011), the WPL (Webb-Pearman-Leuning) correction (Webb et al., 1980), a storage term calculation (Papale et al., 2006; see Appendix A for more details regarding the storage term calculation), spike 30 detection (Papale et al., 2006), gap-filling (marginal distribution sampling method; Reichstein et al., 2005), and nighttime CO2 flux correction (van Gorsel et al., 2009).The details of the gap-filling and partitioning methods are described in the next chapters.
Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.The missing H2O flux (i.e., evapotranspiration, ET) data were gap-filled using the marginal distribution sampling (MDS) method (Reichstein et al., 2005;Hong et al., 2009).This method calculates a median of 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 5 meteorological conditions were 50 W m -2 for the downward shortwave radiation (Rsdn), 2.5°C for the air temperature (Ta), 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 days of 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 less than 56 days, Rsdn was exclusively used following the same approach (i.e., calculating a median 10 of ET under similar Rsdn conditions within a time window).This gap-filling method is used for not only the H2O flux but also the sensible heat and daytime CO2 fluxes.

Modeling of wet canopy evaporation
For estimating the wet canopy evaporation (EWC), 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 15 program.The EWC is estimated as follows: where EWC_Mod is the modeled EWC, σf is the vegetation fraction (i.e., 1-gap fraction); and Ep is the potential evaporation ( VPD / ( 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, cp is the specific heat of air, ga is the aerodynamic 20 conductance, λ is the latent heat of vaporization); Wc is the intercepted canopy water; S is the canopy storage capacity; the exponent n is an empirical coefficient (2/3); and g0 is the architectural conductance (0.5 m s -1 ).
The variable Wc is estimated as: where P is the input total rainfall and D is the drip.When Wc > 0, the canopies are wet.When Wc > S, the drip starts (D > 0). 25 The S is estimated by 0.2 × LAI (leaf area index), and the gap fraction is estimated by exp(-k × LAI), where k varies from 0.3 to 1.5, depending on the species and canopy structure (Jones, 1993, k = 0.75 and 0.485 for the GDK and GCK, respectively.This method only considers the EWC from the canopy by neglecting the EWC from the trunk and stem.Besides, the interception of snow is not considered because the small amount of intercepted snowfall evaporates when the eddy covariance systems 5 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 condition.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.10

Gap-filling and partitioning technique for evapotranspiration: model-stats hybrid method
The currently used MDS is expected to under-and over-estimate 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 (Wc, see Eq. 1) was calculated using the simplified Rutter sparse model; (2) all the 15 missing gaps were filled by the MDS using the data under dry canopy conditions only (i.e., when Wc = 0), which corresponds to the ET under dry canopy condition (ETdry); (3) under wet canopy conditions (i.e., when Wc > 0), the gap-filled data were replaced with the sum of the EWC estimated by the simplified Rutter sparse model (i.e., EWC_Mod) and the ETdry multiplied by 1-(Wc/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 20 transpiration under a dense and closed canopy) and EWC as follows: (1) if the data was gap-filled data (i.e., formerly missing or discarding after quality control), then T was estimated as (1-(Wc/S) n ) ETdry, while the EWC was estimated as EWC_Mod; for the observed ET (ETObs); (2) if the signs of ETObs, ETdry, and EWC_Mod were the same, they were partitioned into T (or EWC) by multiplying ETObs and the ratios of (1-(Wc/S) n ) ETdry (or EWC_Mod) to the sum of (1-(Wc/S) n ) ETdry and EWC_Mod; (3) if not, the T and EWC were estimated by (1-(Wc/S) n ) ETdry and (1-(Wc/S) n ) ETdry, which was then subtracted from ETObs.The procedure 25 regarding the MSH is described in Fig. 1.
[Figure 1  The KoFlux protocol includes three different nighttime corrections methods: the friction velocity (u * ) filtering method (FVF), the light response curve method (LRC), and the modified van Gorsel method (VGF) (van Gorsel et al., 2009;Kang et al., 2014 and2017).These three filtering methods each have their own way of selecting good quality CO2 flux data.The site-specific 5 settings of the individual methods were as follows: 1) the lower u * threshold for the FVF was 0.3 m s -1 for both the GDK and GCK sites based on the dependency of the observed CO2 flux (FCO2_Obs) on the friction velocity at nighttime (for all season, and the higher threshold was not used.Kang et al., 2014), 2) the Michaelis-Menten-type light response curve , where RLRCd is the estimated mean daytime RE, α is the apparent quantum yield, Amax is the canopy scale photosynthetic capacity, and Qt is the total incident shortwave radiation above the canopy; van Gorsel 10 et al., 2009) was used for the LRC, and 3) the peak of the FCO2_Obs that occurred approximately at sunset (Rmax) was directly used for the modified VGF after calculating the median diurnal variation of the CO2 flux for a certain period (Kang et al., 2014 and2017).The modified VGF produces a similar result to that from the original VGF for the sites (see Appendix C for more details regarding the modified VGF).We applied a 30-day moving window to obtain the daily RLRCd and Rmax.
The selected RE data from each filtering method were processed as follows.First, we estimated the parameters in the RE 15 equation (Lloyd-Taylor equation, , where Rref is the reference RE, Tref is the reference temperature (= 10°C), E0 is the activation energy parameter ( o C -1 ), T0 is -46.02°C and Ta is the air temperature ( o C), using the selected observed RE (Lloyd and Taylor, 1994).Second, we replaced the bad (or missing) data with the calculated data using the air temperature and the RE function with the estimated parameters.We estimated Rref using a 30-day moving window which was shifted every 5-days to consider the variations of an RE-controlled by soil moisture and phenology, 20 which is not considered in the Lloyd-Taylor equation.The E0 is constant for each site-year, which is estimated using the generic algorithm proposed by Reichstein et al. (2005) that derives a short-term temperature sensitivity (see Reichstein et al., 2005;Hong et al., 2009 for more details).GPP was calculated by subtracting NEE from RE.

A modification of the moving point test method (MPT)
The objective of the moving point test (MPT) method is to determine the intermediate range of u * where the nighttime CO2 25 fluxes are independent of u * (Gu et al., 2005).The method searches for lower and higher u * thresholds, which are found by statistically testing (i.e., t-test) a group of points with consecutive u * values in a narrow moving window against a reference sample.The original method excludes night data when the median u * were lower than the lower u * threshold, to avoid an underestimation of the CO2 flux due to drainage flow.However, this consideration is inappropriate for hilly terrain sites that are usually affected by drainage flow (e.g., the study sites, GDK and GCK).We modified the original MPT method to apply 30 it to hilly terrain sites by: (1) splitting it into the two time windows, i.e., the time window near sunset (when a drainage has not  2) splitting the time windows, i.e., the first window one and two hours before and after the time of peak occurrence, respectively; and the second window of the time after the first time window.For applying the MPT method, there were two criteria, i.e., the significance level in the t-test (α) and the size (i.e., the number of data) of the moving sample (n).According to Gu et al. (2005), the α is 0.1, and the n is 25.The MPT method was applied for every 10 three months.The details regarding the (modified) MPT are described in the flow charts (Fig. 2) and Gu et al. (2005).
[Figure 2 here] 3 Results 15 3.1 Gap-filling and partitioning results for the H2O flux 3.1.1Validation of the MSH First, we evaluated the latent heat flux under (mostly) wet canopy conditions (λETWC, i.e., λET when Wc/S>2/3) from the model-stats hybrid (MSH) method (λETWC_MSH) against the observed λETWC (λETWC_Obs) at both sites from 2008 to 2010 (Fig. 3).Most of the points are near a one-to-one line.The data scattered away from the one-to-one line are characterized by large 20 aerodynamic conductance (e.g., >100 mm s -1 ) and/or large VPD (e.g., >10 hPa).
The d values for the sites were close to 1 (0.91±0.01 for the GDK and 0.95±0.01for the GCK).Compared to the 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 30 the dry canopy.One of the leading causes of the error in λETWC_MSH was identified as the discrepancy between the time when Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License. the rain occurred, and the tipping bucket was tipped.To validate only EWC, 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 λETWC_MSH and λETWC_Obs show that MSH can provide λETWC reasonably well for the sites.

Comparison between the MDS and the MSH
To evaluate the superiority of the MSH, we filled up the missing λETWC data by using the MDS (λETWC_MDS) and the MSH (λETWC_MSH).Based on the previous research, 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 (ga) and sensible heat 10 advection, Kang et al., 2012).Actually, the average annual MBEs from 2008 to 2010 were -18±6 W m -2 for the GDK site and -15±5 W m -2 for the GCK site, respectively (not shown here).It also should be noted that the λETWC_MSH varied while λETWC_MDS were nearly constant, because (1) the λETWC_Obs rarely existed close to the missing data and (2) the MDS did not consider the effect of ga.
Figure 4 shows the monthly ETs gap-filled by the MDS and MSH methods for the GDK and GCK sites.First, the annual ETs 15 from the MSH method were 16 ~ 41 mm year -1 , which is significantly larger than those from the MDS method, while the random uncertainties in gap-filled annual ETs were approximately 5 mm year -1 for the both sites (quantified according to Finkelstein andSims, 2001, andRichardson andHollinger, 2007).The significant difference was identified in June, July, August and September when it was intensive 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 1,407, 1,323, and 20 1,652 mm in 2008, 2009, and 2010, respectively.Such characteristics are similar to those for the GCK).In addition to taking the missing EWC from the MDS into account, the other advantage of the MSH method is that the observed in ET and by eddy covariance system can be partitioned into transpiration (T) and EWC 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 (ES).Thus, there is more separating the EWC than partitioning the ET.The annual 25 EWC ranged from 53 mm to 82 mm for the GDK and 78 mm to 112 mm for the GCK, which occupies 14 ~ 23% and 14 ~ 19% of the annual ET, respectively.
For quantifying the ES, supplementary eddy covariance (EC), the systems were operated at the floors of the GDK and GCK sites (Kang et al. 2009b).The annual understory ET (~ ES) from 1 June 2008 to 31 May 2009 was 59 mm for the GDK and 43 mm for the GCK, which occupied 16% and 8% of the annual ET, respectively.The decoupling factor (Ω, McNaughton and 30 Jarvis, 1983) at the forest floor was approximately ~ 0.15 for the both sites, which indicates that the ES was controlled primarily by the VPD and surface conductance (gs) rather than Rsdn.This factor also suggests that separating ES from ET using the Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.exponential radiation extinction model to estimate the Rsdn at the forest floor and the relationship between the estimated Rsdn and the ET when the canopy is inactive (Stoy et al., 2006) can be problematic.Considering that the accurate estimation of gs is challenging, a supplementary measurement (e.g., low-level EC, lysimeter, sap-flow measurement, and isotope) is a better approach for estimating ES.Using the ES 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 EWC estimated as 55 mm (15% of ET) and 82 mm (14% of ET), 5 respectively (Fig. 4).
[Figure 4  was statistically the same.The result suggests that there is no drainage and/or the drainage effect is negligible.Meanwhile, 15 during DOY 1-90, those are significantly different, suggesting that the observation for the 2 nd time window cannot represent the actual ecosystem respiration due to the drainage.Many points for the 2 nd time window are near 0 although the air temperature is higher than 0 °C (Fig. 5(b)).The results from the modified MPT are summarized in Table 2.It is hard to find a certain characteristic of the drainage effect, implying that it is a consequence complexly influenced by micrometeorology, phenology, and data availability.Such results require careful examinations from micrometeorological and ecological 20 perspectives using the other independent observations because the modified MPT is a kind of empirical correction.
[Figure 5 here] [Table 1 here]   25 To validate the results from the modified MPT, we compared the results to those from a previous study (Kang et al. 2017).
They developed the site-specific quality control filter to exclude data strongly affected by CO2 advection, which identifies the observations when a high information flow toward the bottom of the slope exists in the dynamical process network based on the observed multi-level CO2 concentrations at the GDK and GCK sites.The filter discards the data (1) when the CO2 drainage is fully generated as the mountain wind prevails (i.e., 21:00 to 9:00, 8:00, and 7:30 for the dormant, transition, and growing 30 seasons, respectively), (2) when u * is lower than the threshold (0.3 m s -1 ) while the drainage flow is under strong development (from 17:00 to 21:00), and (3) until the accumulated CO2 completely dissipates for the downhill (GCK) site (i.e., before noon).
Figure 6 shows the (averaged) annual CO2 budget from the three general nighttime correction methods (see the Chapter 2.4.1) after applying the filter and that from the modified MPT method for the sites.Almost all the results agree with each other within the margin of error, indicating the validity of the modified MPT method.Additionally, we compared the annual net ecosystem production (NEP, = net primary production (NPP) -heterotrophic respiration) with the diameter at breast height (DBH) increment (which is directly related to the NPP) for the GCK site with the single dominant species (i.e., Abies 5 holophylla).The both show similar interannual variability, which also supports the reasonability of the method (Fig. 7).

Applications
In this chapter, we illustrate the advantages of the proposed techniques.Mostly, the benefits are caused by the gap-filling and partitioning of the H2O flux because the model-stats hybrid (MSH) method can take the EWC into account properly and separate them from ET, which has not yet been previously possible.On the other hand, the modified moving point test (MPT) method has an improved applicability because it can be applied regardless of topography.15

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. 8, see Hong et al. (2010) and Grinsted et al. (2004) for more details regarding the wavelet coherence analysis).From one-to the third month period, which was the three-month period in the monsoon season (i.e., the intensive 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 20 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 2006In , 2008In , and 2009, the arrow on the high correlation area pointed left.It means a negative correlation between the two variables, reflecting that the decrease in T caused by the diminishment of Rsdn during the intensive rainy period.In contrast, the arrow pointed right in 2007, indicating a positive correlation.The magnitude of enhanced EWC was greater than that of decreased T at that time and frequency in 2007.Such a positive correlation between ET and rainfall 25 with one-to three-month cycle in 2007 was not reported in the previous study of Hong et al. (2010) 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 EWC in the general gap-filling method).During the monsoon season, the EWC compensates (a portion of) the decreased T, and it can occasionally be balanced (e.g., in 2010).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 CO2 flux to H2O flux), we re-defined the annual ecosystem-level WUE (WUEEco) and the annual canopy-level WUE (WUECanopy) as ΣNEP/ΣET and ΣNPP/ΣT, respectively.For estimating ΣNPP and ΣT simply, we used 0.45 of the ratio of the NPP to GPP for the both sites (Waring et al., 1998), and 0.156 and 0.075 of 10 the ratios of ES to ET for the GDK and GCK, respectively (Kang et al., 2009b).From 2006 to 2010, WUEEco (WUECanopy) ranged from -0.16 (2.17) to 0.32 (2.59) g C•(kg H2O) -1 for the GDK site and from 0.20 (1.93) to 0.38 (2.16) g C•(kg H2O) -1 for the GCK site (Table 3).Considering the increasing trend of NEE and GPP for the GCK site, it can be identified that the interannual variabilities of WUEEco and WUECanopy occurred in opposite directions for the both sites.It was primarily caused by that EWC were enhanced in 2007 and 2010 due to the weakest rainfall intensity and the largest rainfall amount, respectively.15 Overall, such partitioning of the total ET into EWC, T, and ES enables us to understand better how ET responses to environmental changes and the how water cycle is connected to the carbon cycle in a forest ecosystem.

Conclusions 20
There is a common characteristic between the two techniques proposed in this study for gap-filling and partitioning of H2O and CO2 eddy fluxes: two existing methods were merged into a new method.The marginal distribution sampling (MDS) method and the simplified Rutter spars model have merged into the model-stats hybrid (MSH) method and the u * filtering (i.e., moving point test method, MPT) and van Gorsel methods were merged into the modified MPT.Such a strategy strengthens the strength and makes up for the weakness of the original methods.Especially, the modified MPT for nighttime CO2 flux 25 correction substantially improves its applicability, expecting that it will contribute to the standardization of eddy covariance data processing.In this context, such attempt will and must continue.Figure 6: The (averaged) annual CO2 budget (NEE: net ecosystem exchange, GPP: gross primary production, RE: ecosystem respiration) from the three general nighttime correction methods (i.e., u * filtering method, light response curve method, and van Gorsel method) after applying the site-specific filter (adapted from Kang et al., 2017) and that from the modified MPT method for the sites.Error bar indicates the standard deviation of the results from the three general nighttime correction methods.Table A1: Annually integrated net ecosystem exchange (NEE), gross primary production (GPP), and ecosystem respiration (RE), evapotranspiration (ET), and sensible heat flux (H) depending on the estimation of the storage term (profile: estimate of storage term using profile data, single: estimation of storage term using one level (top) data, no storage: not considering storage term; if profile data was not available, we also excluded eddy covariance data for a fair comparison.Accordingly, the values could be slightly different to those in

2. 3
Gap-filling and partitioning methods for the H2O flux 2.3.1 Marginal distribution sampling (MDS) method Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.This parameter is estimated using a plant canopy analyzer (Model LAI-2000; Li-Cor Inc.)) (see Appendix B for more details regarding the parameterizations of S).The input data (and parameters) can be obtained from the flux tower measurement without extra effort.
here] Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.2.4 Correction and partitioning methods for the CO2 flux 2.4.1 General nighttime CO2 flux correction methods Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.yet (completely) manifested) and the time window after the former, (2) applying the MPT method to each time window, (3) comparing the results between the two time windows and determining the existence of CO2 drainage (i.e., the averages of the normalized nighttime CO2 fluxes for the two time windows were significantly different), and (4) excluding data in the second time window if the CO2 drainage is commonly generated and applied and then the u * filtering method using the u * thresholds determined in the previous steps.The best feature of the modified MPT is that the time is split into the two time windows 5 based on van Gorsel et al. (2009): (1) for calculating the median diurnal variation of the CO2 flux and identifying whether the peak of NEE occurred approximately at sunset; ( Figure 5 shows a part of the results from the modified moving point test (MPT) method for the GCK in 2008.During DOY 91-181, the averages of the normalized nighttime CO2 fluxes for the two time windows (i.e., the time window near sunset (when drainage has not yet (completely) manifested (as a 1 st time window) and the time window after the (2 nd time window))

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

Figure 2 :
Figure 2: Flowcharts of the gap-filling and partitioning technique for CO2 flux, the modified moving point test (MPT) method (a) and the original MPT (b; adapted from Gu et al. 2005).

Figure 3 :
Figure 3: Comparison of the latent heat flux under (mostly) wet canopy condition (i.e., Wc/S>2/3 where Wc is the intercepted canopy water and S is the canopy storage capacity) at the GDK (a) and GCK (b) sites: λETWC_Obs indicates the observed latent heat flux under a wet canopy condition (λETWC), while λETWC_MSH indicates the estimated λETWC from the model-stats hybrid method.The dotted line represents the 1:1 line.5

Figure 4 :
Figure 4: Seasonal variation of monthly integrated evapotranspiration (ET) with the gap-filled by the marginal distribution sampling method (ETMDS); the ET gap-filled by the model-stats hybrid (MSH) method (ETMSH), transpiration and wet canopy evaporation partitioned by the MSH method (TMSH and EWC_MSH), for the GDK (a) and GCK (b) sites.ES_Obs indicates soil evaporation measured by the supplementary eddy covariance systems at the floors (adapted from Kang et al. 2009b).5

Figure 5 :
Figure 5: The relationship between the nighttime CO2 flux and air temperature after applying the modified moving point test (MPT) for the GCK in 2008, DOY 91-181 (a) and DOY 1-90 (b).The white color indicates the 1 st time window, whereas black color indicates the 2 nd time window.The solid line indicates the fitting line of Lloyd-Taylor equation using the filtered data of the 1 st time window.

Figure 7 :
Figure 7: The interannual variabilities of the annual net ecosystem production and the annual increment of DBH (diameter at breast height) for the GCK site.

Figure 8 :
Figure 8: Wavelet coherence spectrum of evapotranspiration (ET) with rainfall (P) for the GDK site.A 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 (positive correlation) pointing right, antiphase (negative correlation) pointing left, and P leading ET by 90° pointing straight down).The shaded area indicates the cone of influence where the edge effects might distort the results.5

Figure B1 :
Figure B1: Relationship between canopy storage capacity and plant/leaf area index (the data from TableB1).

Figure C1 :
Figure C1: Seasonal variation of the daily integrated (7-day running mean) net ecosystem exchange (NEE), gross primary production (GPP), and ecosystem respiration (RE) for the GDK and GCK sites in 2008.The VGF and VGFm mean van Gorsel and modified van Gorsel methods, respectively.The site-specific quality control filter for eliminating the drainage-affected data (Kang et al., 2017) was not applied.The values in parenthesizes mean the annually integrated NEE, GPP, and RE (unit: g C m -2 y -1 ).5

Table 1 :
Statistical parameters for 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 2 are from the linear regression analysis.Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.

Table 2 :
Summary results generated from the modified moving point test (MPT) for the GDK and GCK sites.u * L and u * H indicate lower and higher u * thresholds, 1 st and 2 nd indicate the 1 st and 2 nd time windows, respectively.Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.

Table B1 :
Review of the parameters for wet canopy evaporation modeling in the previous studies.Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-247Manuscript under review for journal Biogeosciences Discussion started: 14 June 2017 c Author(s) 2017.CC BY 4.0 License.