Spatial variations in evapotranspiration over East Asian forest sites. II. Surface conductance and aerodynamic conductance

Evapotranspiration over forest surfaces is mainly constrained by environmental and forest structural components through their influence on surface conductance (Gs) and aerodynamic conductance (Ga). Tower based eddy covariance data from 16 forest sites in East Asia (2°S to 64°N) were used to examine the sensitivity of Matsumoto and Nakai models for predicting Gs and Ga, respectively. Daytime dry-canopy data for the growing season were used in this study. Comparisons between model predictions and observed Gs and Ga showed good agreement, suggesting that the models were suitable for predicting Gs and Ga with reasonable accuracy. However, the model for Gs was unable to predict Gs accurately when soil water content was low (∼10%). In this circumstance, effective soil water content and a more comprehensive method for modelling the soil water content function must be used. Gs in East Asia was largely depended on vapour pressure deficit and secondarily on soil water content. Ga was largely affected by leaf area index compared with stand density.


Abstract:
Evapotranspiration over forest surfaces is mainly constrained by environmental and forest structural components through their influence on surface conductance (G s ) and aerodynamic conductance (G a ). Tower based eddy covariance data from 16 forest sites in East Asia (2°S to 64°N) were used to examine the sensitivity of Matsumoto and Nakai models for predicting G s and G a , respectively. Daytime dry-canopy data for the growing season were used in this study. Comparisons between model predictions and observed G s and G a showed good agreement, suggesting that the models were suitable for predicting G s and G a with reasonable accuracy. However, the model for G s was unable to predict G s accurately when soil water content was low (~10%). In this circumstance, effective soil water content and a more comprehensive method for modelling the soil water content function must be used. G s in East Asia was largely depended on vapour pressure deficit and secondarily on soil water content. G a was largely affected by leaf area index compared with stand density.

INTRODUCTION
Evapotranspiration (ET) is a major component influencing the global hydrologic process and thus climate. Surface conductance, G s , and aerodynamic conductance, G a , are important variables for determining the ET over vegetated surfaces Khatun et al., 2011). G s is known to be controlled by several environmental components (e.g., radiation, temperature, vapour pressure deficit, and soil water content) . G a is governed by wind speed and surface aerodynamic properties, such as roughness length (z 0 ) and zero-plane displacement (d), which are the functions of forest structural components (e.g., stand density, ρ s , and leaf area index, LAI) . It is therefore important to quantify how different environmental and forest structural components control G s and G a to understand the relationship of ET to these components. G s is usually calculated using the inverse form of the Penman-Monteith equation (Monteith, 1965), and G a has generally been computed from wind speed and friction velocity based on observed flux data. However, the accurate calculation of G s and G a using flux data is sometimes difficult because of the heterogeneity of the landscape and the large number of controlling factors involved, including climate, plant biophysics, soil properties and topography.
A number of empirical models have been proposed to calculate G s for many stands of various vegetation categories, including grasslands, crops and broadleaved and coniferous forests (e.g., Jarvis, 1976;Ball et al., 1987;Leuning, 1995). On the basis of the Jarvis model, Matsumoto et al. (2008) developed a pooled model of three temperate and two sub-arctic forests that allows the expression of G s in different forests under different climatic zones using a common parameter set.
Many different approaches have been developed to calculate d and z 0 (e.g., Pinard, 2000;Raupach, 1994;Choudhury and Monteith, 1988). However, existing approaches generally utilise highly simplified representations of vegetation structure. Nakai et al. (2008) proposed a new model to calculate d and z 0 using forest structure data such as LAI (m 2 m −2 ) and ρ s (trees ha −1 ).
However, the Matsumoto and Nakai models that predict G s and G a are derived from small numbers of measurements. In this study, we compiled the measurement data from 16 forest sites in East Asia to evaluate the effectiveness of the Matsumoto and Nakai models in relation to our data sets. We also examined the dependence of G s and G a on environmental and forest structural components using these Correspondence to: Rehana Khatun, Lab. of Forest Meteorology and Hydrology, Graduate School of Bioagricultural Sciences, Nagoya University, Japan. E-mail: rehanakhatun_2008@yahoo.co.jp ©2011, Japan Society of Hydrology and Water Resources. models. This information will be useful for land surface parameterisation in regional climate simulations, particularly in East Asia.

MATERIALS AND METHODS
Site description, measurement system, quality control and data selection for turbulent fluxes Data were collected from 16 forest sites in East Asia that were distributed geographically from 2°S to 64°N latitude and 98°E to 142°E longitude (see Figure S1 for the geographic and climatic positions in detail). The study included one tropical rain forest (PDF), four tropical monsoon forests (SKR, MKL, MMP and KMW), one subtropical monsoon forest (QYZ), two warm-temperate forests (SMF and GDK), one temperate forest influenced by the monsoon (CBS), four cool-temperate forests (TMK, MBF, MMF and TSE), and three sub-arctic forests (SKT, YLF and TUR).
An eddy covariance and a meteorological measurement system were used at each site to measure sensible heat flux (H), latent heat flux (λE), friction velocity (u * ) and basic environmental components. The principal investigators performed quality controls at each site. We collected 30and 60-min averages of flux and meteorological data for this study. Data for daytime (net radiation >0 W m −2 ) and dry-canopy conditions (excluding data collected during rainfall events and within 10 h after these events) for the growing season were used for the analysis. After excluding the wet-canopy data, more than 60% of daytime data were available for most of the measurement years. To avoid variations in the closure of energy balance among the sites, closure was forced to 1 by the Bowen ratio closure method.
The site description (Table S1), eddy covariance instruments (Table S2) and quality control and data selection for turbulent fluxes (Text S1) are described with more details in the supplements. The calculation of G s and G a based on observed data is described in Text S2.

Estimation of G s and G a
According to the Jarvis-type model (Jarvis, 1976), G s is a function of several environmental variables and is expressed as where G smax is the maximum surface conductance and f 1 (Q), f 2 (D), f 3 (T), and f 4 (θ) are the functions of the photosynthetic photon flux density (Q, μmol m −2 s −1 ), vapour pressure deficit (D, kPa), air temperature (T, °C) and volumetric soil water content (θ, %), respectively. Each function represents the influence of that factor on G s and ranged from 0 to 1.
According to Matsumoto et al. (2008), the functions are expressed as follows: (2) where D 0.5 , T min , T opt , T max , θ min and k 1 -k 4 are the fitting parameters; Q max and θ max are fixed parameters. T min , T opt and T max are the minimum, optimum and maximum T values, respectively; θ min and θ max are the minimum and maximum θ values; and Q max is the maximum Q value. Detailed descriptions of all fitting and fixed parameters were provided by Matsumoto et al. (2008). According to the pooled model, the values of the fitting parameters are k 1 = 310, k 2 = 3.0, k 3 = 0.10, k 4 = 20, D 0.5 = 1.25, T min = 2, T opt = 22, T max = 37 and θ min = 10 and fixed parameters are Q max = 2100 and θ max = 50; G smax = 0.037 m s −1 , when LAI ≥ 3; G smax = G sMAX (LAI/6 + 0.5) when 0 < LAI < 3; and G sMAX is the maximum G s with highest LAI.
The effects of each variable on the variation of G s were calculated as the contribution of each function (α x = β x − ω a ) to modelling accuracy using the root mean-square error (RMSE) value (Matsumoto et al., 2005). α x is the contribution index (m s −1 ) of one function, i.e., f(x) to the variability of G s , β x is the RMSE between the observed and estimated values in which a function for a certain variable was excluded, and ω a is the RMSE between the observed and estimated values containing all of the functions in Equation (1). More description is in Text S3.
G a was calculated as (Businger, 1956) where z is the measurement height (m), U z is wind speed at height z (m s −1 ), and k is von Karman's constant (0.40). According to the Nakai et al. (2008) model, d and z 0 are expressed as (7)   (8) where d/h and z 0 /h are the zero-plane displacement and roughness length, respectively, normalised by canopy height h, and α and β are the fitting parameters.

Surface conductance
pooled model clearly expresses the response characteristics of stomata to each variable as of the study of Matsumoto et al. (2005). However, at two sites (SKR and SKT) with low θ, G sobs exceeded the potential curve represented by f(θ). This suggests that f(θ) could not accurately represent the G sobs for low-water conditions in the shallow soil layer, whereas the other functions of this model are potentially applicable to our data set.
The linear regression between modelled G s (G smod ) calculated by the pooled model and G sobs (Figure 2a) gave a slope of 1.06 with a high level of significance (P < 0.005) and RMSE of 0.005 m s −1 . This result indicates that the pooled model is capable of predicting G s . However, the relatively low R 2 value (0.51) indicates that G smod did not always fit the calculated value well. This result might be explained by the large underestimation of G s by the pooled model at tropical monsoon site SKR and boreal site SKT and overestimation of G s at cool-temperate sites MMF and TSE.
The poor agreement at SKR and SKT might be partially attributable to the use of top layer θ in calculating the f(θ). The low water-holding capacity of the sandy soil at SKR and SKT may cause a low-water condition in the top soil layer that greatly affects model values of G s . However, G sobs at these sites was not affected by low θ (Figure 2). Trees might be able to use deeper soil water when the top layer becomes dry (e.g., Tanaka et al., 2004;Li et al., 2006). On the other hand, the pooled model overestimated G s at two cool-temperate sites (MMF and TSE). At these sites, the G smod values were approximately double than those of the G sobs . Higher θ values at these sites resulted in higher G s model output. However, in real situation, θ had lower impact on the actual G s . These sites are belonging to maritime climate and are characterized by low temperature and high humidity resulting in low atmospheric evaporative demand.
The lower values of latent heat flux and higher values of Bowen ratio due to low atmospheric evaporative demand might affect the actual G s (see Equation (1) in Text S2). Without these exceptional four sites, the linear regression between G smod and G sobs gave a slope of 0.90 and R 2 of 0.75 with P < 0.0005 and RMSE of 0.003 m s −1 . This high correlation illustrates that in absence of the water stress (θ 10%) condition or high water content (θ >40%) condition with low evaporative demand, the pooled model is suitable for the prediction of G s . Figure 2b shows the contribution (α x ) of each function to the variability of G s . Among all of the variables, D was the dominant component of variability in G s , and the order of α x for each function was The values of f(D) were reduced by the higher values of D at most of the sites and were a greater controlling factor for G s . The exceptions were found at cool-temperate sites, where the values of f(D) were nearly 1 and not a controlling factor for G s . At tropical, warm-temperate and some boreal sites, where θ <30%, f(θ) was the greater controlling factor. The values of f(Q) were >0.80 at all sites except two boreal sites. At boreal sites, G s was quite affected by Q. At most of the studied sites, T was between 17°C and 27°C [f(T) >0.90] and did not act as a controlling factor. However, at some cool-temperate and boreal sites, T was between 12°C and 15°C and exerted a control on G s . Therefore, G s in East Asian forests was not directly determined by Q and T. Even if the functions of Q and T were excluded from the pooled model, R 2 and P values became slightly larger than those of the complete model. This was probably the result of the  high correlation between T and D (R 2 = 0.74), indicating that the effect of T was incorporated into D and the narrow range of growing season mean Q and T. However, without the functions of Q and T, RMSE was increased, indicating that estimation of G s without these functions was less satisfactory at our sites. Therefore, including all functions in the model enhances the applicability of the model at studied sites. Because of previous studies on the response of G s to environmental components (e.g., Ogink-Hendriks, 1995;Matsumoto et al., 2005) have been limited in terms of latitude ranges, the results of some single-site studies on the response of G s to environmental factors based on seasonal variation have discussed here. In those studies, radiation was the dominant factor for the seasonal variation in G s , with results based on hourly values (Ogink-Hendriks, 1995) and daily mean values (Matsumoto et al., 2005) for the growing season. Both hourly and daily mean values can exhibit low radiation that limits G s . In our study, the growing season mean values of Q were mostly >700 μmol m −2 s −1 [f(Q) >0.80] and were not a controlling factor. However, the higher T in the growing season accelerated D, which might have greatly determined dry-canopy G s .

Aerodynamic conductance
Normalised displacement height (d/h) (Equation 7), and roughness length (z 0 /h) (Equation 8) were calculated for 10 forest sites (Table S1) because of the lack of ρ s and LAI data at other sites. The values of d and z 0 estimated by Nakai-model (d mod and z 0mod ) showed a good agreement with those of observed d and z 0 (d obs and z 0obs ) ( Figure S2 and detailed in Text S4). However, the validity test of Nakaimodel for estimating d and z 0 was limited at six sites only (two tropical, two temperate and two boreal sites) because of the lack of wind profile data at other sites.
G a estimated by Nakai-model (G amod ) was compared with observed G a (G aobs ) calculated by measured friction velocity and wind speed. Linear regression between G amod and G aobs (Figure 3a) gave a slope of 0.91 and R 2 of 0.70 with a high level of significance (P < 0.002), and RMSE of 0.015 m s −1 . This result indicates that the Nakai-model is highly satisfactory for the prediction of d and z 0 , and is thus applicable to the prediction of G a .
To quantify the contribution (α x ) of ρ s and LAI to the variability in G a , RMSE was calculated between G aobs and G amod using the estimated values of d/h and z 0 /h, where ρ s and LAI were excluded from Equation 7 one by one (Text S5). As shown in Figure 3b, LAI made a larger contribution than did ρ s to the variation in G a , because of the close relationship of the roughness parameters (d/h and z 0 /h) to LAI compared with ρ s , i.e., d/h increased and z 0 /h decreased clearly with the increase in LAI, however, d/h increased and z 0 /h decreased with the increase in ρ s upto certain limit and after that changing only slightly as ρ s increased ( Figure S3).

SUMMARY AND CONCLUSIONS
Our results indicate that both the Matsumoto and Nakai models can be applied to broad climate and forest types using growing season mean data. However, G s was poorly modelled at extremely low soil water condition (θ of approximately 10%). The use of appropriate θ values and a more comprehensive method for modelling the soil water function are needed under this condition. G s depended on the environmental components in the order D > θ > Q > T. The dependence of G a on LAI was higher than that on ρ s . This study is the first to predict G s and G a using the Matsumoto and Nakai models in East Asian forests. The results provide useful information for predicting G s and G a and indicate that the estimation of these parameters is possible in the absence of flux data at a larger scale.

SUPPLEMENTS
Supplement 1. This includes: Table S1. Description of the study sites Table S2. Measurement system, height, year, closure rate and growing season length Text S1. Description of quality control and data selection for turbulent fluxes Text S2. Calculation of G a and G s based on observed data Text S3. Calculation of dependence of G s on environmental components Text S4. Calculation of zero-plane displacement (d) and roughness length (z 0 ) Text S5. Calculation of dependence of G a on forest structural components Figure S1. Locations of the study sites: (a) geographic and (b) climatic (based on the annual sum precipitation and annual mean temperature). Figure S2. Relationship between: (a) estimated (d mod ) and observed (d obs ) zero-plane displacement (d) and (b) estimated (z 0mod ) and observed (z 0obs ) roughness length (z 0 ). Figure S3. Relationship between normalized zero-plane displacement (d/h), roughness length (z 0 /h) and (a) stand density (ρ s ), (b) leaf area index (LAI).