**Development of Hybrid Method for the Modeling of Evaporation and Evapotranspiration**

Sungwon Kim *Dongyang University, Republic of Korea* 

## **1. Introduction**

350 Evapotranspiration – Remote Sensing and Modeling

Tasumi, M.; Trezza, R.; Allen, R.G. & Wright, J.L. (2005). Operational Aspects of Satellite

Taylor, S.A. & Ashcroft, G.L. (1972). *Physical Edaphology. The Physics of Irrigated and Non-*

Thornthwaite, C.W. (1948). An approach toward a rational classification of climate.

Turc, L. (1961). Evaluation de Besoins en Eau d'Irrigation, ET Potentielle. *Ann. Agron*., Vol.

Van Genuchten, M.Th. (1980). A Closed-Form Equation for Predicting the Hydraulic

Van Walsum, P.E.V.; Veldhuizen, A.A.; van Bakel, P.J.T.; van der Bolt, F.J.E.; Dik, P.E.;

Wang, X.; Xie, H.; Guan, H. & Zhou, X. (2003). Different Responses of MODIS-Derived

Wapener, T. (2003) Evaluation of Catchment Models. *Hydrological Processes*, Vol. 17, pp.

Xu, Y. & Beekman, H.E. (Eds.) (2003). *Groundwater Recharge Estimation in Southern Africa*,

Yang, D.; Zhang, T.; Zhang, K.; Greenwood, D.; Hammond, J.P. & White P.J. (2009). An

Zhang, L.; Dawes, W.R. & Hatton, T.J. (1996). Modelling Hydrologic Processes Using a

Zhang, L.; Dawes, W.R. & Walker, G.R. (2001). The Response of Mean Annual

Zhou, M.C.; Ishidaira, H.; Hapuarachchi, H.P.; Magome, J.; Kiem, A.S. & Takeuchi, K. (2006).

Mokeng River Basin. *Journal of Hydrology*, Vol. 327, pp. 151-173.

ISBN 92-9220-000-3, UNESCO IHP Series no. 64, UNESCO, Paris.

MOBILHY. *Journal of Hydrology*, Vol. 185, pp. 147-169.

*Environmental Modelling & Software*, Vol. 24, pp. 1209-1222.

Conductivity of Unsaturated Soils. *Soil Science Society of America Journal*, Vol. 44(5),

Groenendijk, P.; Querner, E.P. & Smit, M.R.F. (2004). *SIMGRO 5.0.1, Theory and Model Implementation,* Report No. 913.1, Alterra, Wageningen, The Netherlands. Viviroli, D.; Zappa, M.; Gurtz, J. & Weingartner, R. (2009). An Introduction to the

Hydrological Modelling System PREVAH and its Pre- and Post-Processing Tools.

NDVI to Root-Zone Soil Moisture in Semi-Arid and Humid Regions. *Journal of* 

Easily Implemented Agro-hydrological Procedure with Dynamic Root Simulation for Water Transfer in the Crop–Soil System: Validation and Application. *Journal of* 

Biophysically Based Model – Application of WAVES to FIFE and HAPEX-

Evapotranspiration to Vegetation Changes at Catchment Scale. *Water Resources* 

Estimating Potential Evapoutranspiration Using Shuttleworth-Wallace Model and NOAA-AVHRR NDVI Data to Feed a Distributed Hydrological Model over the

*and Drainage Systems*, Vol. 19, pp. 355-376.

*Geograph. Rev*., Vol. 38, pp. 55.

*Hydrology*, Vol. 340, pp. 12-24.

*Hydrology*, Vol. 370, 177–190

*Research,* Vol. 37, pp. 701-708.

12, pp. 13-49.

pp. 892–898.

3375-3378.

*Irrigated Soils*. Freeman, San Francisco, pp. 533.

Based Energy Balance Models for Irrigated Crops in the Semi-Arid US. *Irrigation* 

Evaporation is the process whereby liquid water is converted to water vapor and removed from the evaporating surface. In hydrological practice, the estimation of evaporation can be achieved by direct or indirect methods. Direct method is based on the field measurements. Evaporation pan have also been used and compared to estimate evaporation by researchers (Choudhury, 1999; McKenzie and Craig, 2001; Vallet-Coulomb et al., 2001). The Class A evaporation pan is one of the most widely used instruments for the measurements of evaporation from a free water surface. The pan evaporation (PE) is widely used to estimate the evaporation of lakes and reservoirs (Finch, 2001). Many researchers have tried to estimate the evaporation through the indirect methods using the climatic variables, but some of these techniques require the data which cannot be easily obtained (Rosenberry et al., 2007).

Evapotranspiration (ET) is the sum of volume of water used by vegetation, evaporated from the soil, and intercepted precipitation (Singh, 1988). ET plays an important role in our environment at global, regional, and local scales. ET is observed using a lysimeter directly or can be calculated using the water balance method or the climatic variables indirectly. Because the measurements of ET using a lysimeter directly, however, requires much unnecessary time and needs correct and careful experience, it is not always possible in field measurements. Thus, an empirical approach based on the climatic variables is generally used to calculate ET (Penman, 1948; Allen et al*.*, 1989). In the early 1970s, the Food and Agricultural Organization of the United Nations (FAO), Rome, developed practical procedures to calculate the crop water requirements (Doorenbos & Pruitt, 1977), which have become the widely accepted standard for irrigation studies. A common practice for estimating ET from a well-watered agricultural crop is to calculate the reference crop ET such as the grass reference ET (ETo) or the alfalfa reference ET (ETr) from a standard surface and to apply an appropriate empirical crop coefficient, which accounts for the difference between the standard surface and the crop ET.

Recently, the outstanding results using the neural networks model in the fields of PE and ET modeling have been obtained (Bruton et al., 2000; Sudheer et al., 2003; Trajkovic et al., 2003; Trajkovic, 2005; Keskin and Terzi, 2006; Kisi, 2006; Kisi, 2007; Zanetti et al., 2007; Jain et al., 2008; Kumar et al., 2008; Landeras et al., 2008; Kisi, 2009; Kumar et al., 2009; Tabari et al., 2009; Chang et al., 2010; Guven and Kisi, 2011). Sudheer et al*.* (2002) investigated the prediction of Class A PE using the neural networks model. They used the neural networks model for the evaporation process using proper combinations of the observed climate variables such as temperature, relative humidity, sunshine duration, and wind speed for the neural networks model. Shiri and Kisi (2011) investigated the ability of genetic programming (GP) to improve the accuracy of daily evaporation estimation. They used proper combinations of air temperature, sunshine hours, wind speed, relative humidity, and solar radiation for GP. Kumar et al. (2002) developed the neural networks models to calculate the daily ET. They used proper combinations of the observed climatic variables such as solar radiation, temperature, relative humidity, and wind speed for the neural networks models. Kisi & Ozturk (2007) used the neuro-fuzzy models to calculate FAO-56 PM ETo using the observed climatic variables. They used proper combinations of the observed climatic variables such as air temperature, solar radiation, wind speed, and relative humidity for the neuro-fuzzy models. Kim & Kim (2008) developed the neural networks model embedding the genetic algorithm for the modeling of the daily PE and ETr simultaneously, and constructed the optimal neural networks model using the uncertainty analysis of the input layer nodes/variables. Furthermore, they suggested the 2-dimensional and 3-dimensional maps for PE and ETr to provide the reference data for irrigation and drainage system, Republic of Korea. And, the recent researches combining the stochastic models and the neural networks models in the fields of hydrology and water resources have been accomplished. Mishra et al. (2007) developed a hybrid model, which combined a linear stochastic model and a nonlinear neural networks model, for drought forecasting. Kim (2011) investigated the modeling of the monthly PE and ETr simultaneously using the specific method, which combined the stochastic model with the neural networks models.

The purpose of this study is to develop the hybrid method for the modeling of the monthly PE and FAO-56 PM ETo simultaneously. The hybrid method represents the combination of Univariate Seasonal periodic autoregressive moving average (PARMA) model and support vector machine neural networks model (SVM-NNM). For this research, first, the stochastic model, Univariate Seasonal PARMA(1,1) model, is used for the generation of the reliable data, which are considered as the training dataset. Therefore, the observed data are considered as the testing dataset. Second, the neural networks model, SVM-NNM, is used for the modeling of the monthly PE and FAO-56 PM ETo simultaneously. Homogeneity evaluation using the One-way ANOVA and Mann-Whitney *U* test, furthermore, is carried out for the observed and calculated PE and FAO-56 PM ETo data. And, the correlation relationship between the observed PE and FAO-56 PM ETo data can be derived using the bivariate linear regression analysis model (BLRAM), respectively.

### **2. Calculation of FAO-56 PM ETo**

Penman (1948) combination method links evaporation dynamics with the flux of net radiation and aerodynamic transport characteristics of the natural surface. Based on the observations that latent heat transfer in plant stem is influenced not only by these abiotic factors, Monteith (1965) introduced a surface conductance term that accounted for the response of leaf stomata to its hydrologic environment. This modified form of the Penman-Monteith (PM) ETo model. Jensen et al. (1990) measured ETo using the lysimeters at 11 stations located in the different climatic zones of various regions around the world. They

2009; Chang et al., 2010; Guven and Kisi, 2011). Sudheer et al*.* (2002) investigated the prediction of Class A PE using the neural networks model. They used the neural networks model for the evaporation process using proper combinations of the observed climate variables such as temperature, relative humidity, sunshine duration, and wind speed for the neural networks model. Shiri and Kisi (2011) investigated the ability of genetic programming (GP) to improve the accuracy of daily evaporation estimation. They used proper combinations of air temperature, sunshine hours, wind speed, relative humidity, and solar radiation for GP. Kumar et al. (2002) developed the neural networks models to calculate the daily ET. They used proper combinations of the observed climatic variables such as solar radiation, temperature, relative humidity, and wind speed for the neural networks models. Kisi & Ozturk (2007) used the neuro-fuzzy models to calculate FAO-56 PM ETo using the observed climatic variables. They used proper combinations of the observed climatic variables such as air temperature, solar radiation, wind speed, and relative humidity for the neuro-fuzzy models. Kim & Kim (2008) developed the neural networks model embedding the genetic algorithm for the modeling of the daily PE and ETr simultaneously, and constructed the optimal neural networks model using the uncertainty analysis of the input layer nodes/variables. Furthermore, they suggested the 2-dimensional and 3-dimensional maps for PE and ETr to provide the reference data for irrigation and drainage system, Republic of Korea. And, the recent researches combining the stochastic models and the neural networks models in the fields of hydrology and water resources have been accomplished. Mishra et al. (2007) developed a hybrid model, which combined a linear stochastic model and a nonlinear neural networks model, for drought forecasting. Kim (2011) investigated the modeling of the monthly PE and ETr simultaneously using the specific method, which combined the stochastic model with the

The purpose of this study is to develop the hybrid method for the modeling of the monthly PE and FAO-56 PM ETo simultaneously. The hybrid method represents the combination of Univariate Seasonal periodic autoregressive moving average (PARMA) model and support vector machine neural networks model (SVM-NNM). For this research, first, the stochastic model, Univariate Seasonal PARMA(1,1) model, is used for the generation of the reliable data, which are considered as the training dataset. Therefore, the observed data are considered as the testing dataset. Second, the neural networks model, SVM-NNM, is used for the modeling of the monthly PE and FAO-56 PM ETo simultaneously. Homogeneity evaluation using the One-way ANOVA and Mann-Whitney *U* test, furthermore, is carried out for the observed and calculated PE and FAO-56 PM ETo data. And, the correlation relationship between the observed PE and FAO-56 PM ETo data can be derived using the

Penman (1948) combination method links evaporation dynamics with the flux of net radiation and aerodynamic transport characteristics of the natural surface. Based on the observations that latent heat transfer in plant stem is influenced not only by these abiotic factors, Monteith (1965) introduced a surface conductance term that accounted for the response of leaf stomata to its hydrologic environment. This modified form of the Penman-Monteith (PM) ETo model. Jensen et al. (1990) measured ETo using the lysimeters at 11 stations located in the different climatic zones of various regions around the world. They

bivariate linear regression analysis model (BLRAM), respectively.

neural networks models.

**2. Calculation of FAO-56 PM ETo** 

compared the results of the lysimeters with those of 20 different empirical equations and methodologies for ETo measurements. It was found that PM ETo model showed the optimal results over all the climatic zones. If the observed/measured data for ETo does not exist, therefore, PM ETo model can be considered as a standard methodology to calculate ETo. In Gwangju and Haenam stations which were selected for this study, there are no observed data for ETo using a lysimeter. The data calculated using PM ETo model can be assumed as the observed ETo, whose reliability was verified by many previous studies. All calculation procedures as used in PM ETo model are based on the FAO guidelines as laid down in the publication No. 56 of the Irrigation and Drainage Series of FAO ″Crop Evapotranspiration– Guidelines for Computing Crop Water Requirements″ (1998). Therefore, FAO-56 PM ETo equation means PM ETo equation suggested by the Irrigation and Drainage Paper No. 56, FAO. FAO-56 PM ETo equation is given by Allen et al. (1998) and can be shown as the following equation (1).

$$\text{FAO-56 PM ET}\_{\text{o}} = \frac{0.408\Delta (\text{R}\_{\text{n}} - \text{G}) + \chi (900/(\text{T} + 273)) \text{u}\_{2} (\text{e}\_{\text{s}} - \text{e}\_{\text{a}})}{\Delta + \chi (1 + 0.34 \text{u}\_{2})} \tag{1}$$

where FAO-56 PM ETo = the grass reference evapotranspiration (mm/day); Rn = the net radiation at the crop surface (MJ/m2 day); G = the soil heat flux density (MJ/m2 day); T = the mean air temperature at 2m height (◦C); u2 = the wind speed at 2m height (m/sec); es = the saturation vapor pressure (kPa); ea = the actual vapor pressure (kPa); es - ea = the saturation vapor pressure deficit (kPa); Δ = the slope vapor pressure curve (kPa/◦C); and γ = the psychometric constant (kPa/◦C). FAO CROPWAT 8.0 computer program has been used to calculate FAO-56 PM ETo and extraterrestrial radiation (Ra). FAO CROPWAT 8.0 computer program allows the user to enter the climatic data available including maximum temperature (Tmax), minimum temperature (Tmin), mean relative humidity (RHmean), mean wind speed (WSmean), and sunshine duration (SD) for calculating FAO-56 PM ETo. On the base of climatic data available, FAO CROPWAT 8.0 computer program calculates the solar radiation reaching soil surface. Fig. 1(a)-(b) show the calculation of FAO-56 PM ETo using FAO CROPWAT 8.0 computer program in Gwangju and Haenam stations, respectively.

Fig. 1. Calculation of FAO-56 PM ETo using FAO CROPWAT 8.0 Computer Program

## **3. Stochastic model**

#### **3.1 Univariate seasonal periodic autoregressive moving average model**

Stationary ARMA models have been widely applied in stochastic hydrology for modeling of annual time series where the mean, variance, and the correlation structure do not depend on time. For seasonal hydrologic time series, such as monthly series, seasonal statistics including the mean and standard deviation may be reproduced by a stationary ARMA model by means of standardizing the underlying seasonal series. Hydrologic time series such as monthly streamflows, PE, and FAO-56 PM ETo are usually characterized by different dependence structure depending on the season (Salas, 1993). One may extend Univariate Seasonal periodic autoregressive (PAR) model to include periodic moving average (MA) parameters. Such a model is Univariate Seasonal periodic autoregressive moving average (PARMA) model and is expressed as Univariate Seasonal PARMA(p,q) model. The stochastic models are generally simple to use. When the errors, however, happen in model identification and parameter estimation, the generated data using the stochastic models cannot reconstruct the statistical properties of the observed data exactly. Furthermore, the high-order PARMA(p,q) models have generally many parameters, and the calculation process is much complex (Salas et al., 1980). In this study, the author determined in advance 4 kinds of Univariate Seasonal PARMA(p,q) models including PARMA(1,1), PARMA(1,2), PARMA (2,1), and PARMA(2,2), which are the low-order models and contain the seasonal properties. In general, the low-order Univariate Seasonal PARMA(p,q) models are useful for the periodic hydrologic time series modeling (Salas et al., 1982). Furthermore, the author generated 100 years data in advance using each Univariate Seasonal PARMA(p,q) model for the climatic variables of the neural networks models, respectively. As a result, the author selected Univariate Seasonal PARMA(1,1) model, which shows the best statistical properties and is simple in parameter estimation. Univariate Seasonal PARMA(1,1) model has been applied to monthly streamflow time series from the previous studies (Tao and Delleur, 1976; Hirsch, 1979), and is shown as the following equation (2).

$$\mathbf{v}\_{\mathbf{v},\tau} = \boldsymbol{\mu}\_{\tau} + \boldsymbol{\Phi}\_{\mathbf{1},\tau} (\mathbf{y}\_{\mathbf{v},\tau-1} - \boldsymbol{\mu}\_{\tau-1}) + \boldsymbol{\varepsilon}\_{\mathbf{v},\tau} - \boldsymbol{\theta}\_{\mathbf{1},\tau} \boldsymbol{\varepsilon}\_{\mathbf{v},\tau-1} \tag{2}$$

where yv, / yv, 1 = the monthly PE and FAO-56 PM ETo for year= v and month= / 1 ; / <sup>1</sup> = the means for month= / 1 ; 1, = the seasonal autoregressive parameter for month= ; 1, = the seasonal moving average parameter for month= ; v, v, 1 / = uncorrelated noise terms; v = year; = month (1,2,…, ); and =12. Furthermore, Univariate Seasonal PARMA(2,1), PARMA(2,2) models, and more complex multiplicative PARMA(p,q) models may be needed for hydrologic modeling and simulation when the preservation of both the seasonal and the annual statistics is desired (Salas and Abdelmohsen, 1993).

#### **3.2 Construction of Univariate Seasonal PARMA(p,q) model**

The author used Univariate Seasonal PARMA(1,1) model to generate the sufficient training dataset, and obtained two generated samples. They included the input nodes/variables including mean temperature (Tmean), maximum temperature (Tmax), minimum temperature (Tmin), mean dew point temperature (DPmean), minimum relative humidity (RHmin), mean relative humidity (RHmean), mean wind speed (WSmean), maximum wind speed (WSmax), and sunshine duration (SD) in mean values and the output nodes/variables including PE and

Stationary ARMA models have been widely applied in stochastic hydrology for modeling of annual time series where the mean, variance, and the correlation structure do not depend on time. For seasonal hydrologic time series, such as monthly series, seasonal statistics including the mean and standard deviation may be reproduced by a stationary ARMA model by means of standardizing the underlying seasonal series. Hydrologic time series such as monthly streamflows, PE, and FAO-56 PM ETo are usually characterized by different dependence structure depending on the season (Salas, 1993). One may extend Univariate Seasonal periodic autoregressive (PAR) model to include periodic moving average (MA) parameters. Such a model is Univariate Seasonal periodic autoregressive moving average (PARMA) model and is expressed as Univariate Seasonal PARMA(p,q) model. The stochastic models are generally simple to use. When the errors, however, happen in model identification and parameter estimation, the generated data using the stochastic models cannot reconstruct the statistical properties of the observed data exactly. Furthermore, the high-order PARMA(p,q) models have generally many parameters, and the calculation process is much complex (Salas et al., 1980). In this study, the author determined in advance 4 kinds of Univariate Seasonal PARMA(p,q) models including PARMA(1,1), PARMA(1,2), PARMA (2,1), and PARMA(2,2), which are the low-order models and contain the seasonal properties. In general, the low-order Univariate Seasonal PARMA(p,q) models are useful for the periodic hydrologic time series modeling (Salas et al., 1982). Furthermore, the author generated 100 years data in advance using each Univariate Seasonal PARMA(p,q) model for the climatic variables of the neural networks models, respectively. As a result, the author selected Univariate Seasonal PARMA(1,1) model, which shows the best statistical properties and is simple in parameter estimation. Univariate Seasonal PARMA(1,1) model has been applied to monthly streamflow time series from the previous studies (Tao and Delleur, 1976; Hirsch, 1979), and is shown as the following equation (2).

y (y ) v, 1, v, 1 1 v, 1, v, 1 (2)

where yv, / yv, 1 = the monthly PE and FAO-56 PM ETo for year= v and month= / 1 ; / <sup>1</sup> = the means for month= / 1 ; 1, = the seasonal autoregressive parameter for month= ; 1, = the seasonal moving average parameter for month= ; v, v, 1 / = uncorrelated noise terms; v = year; = month (1,2,…, ); and =12. Furthermore, Univariate Seasonal PARMA(2,1), PARMA(2,2) models, and more complex multiplicative PARMA(p,q) models may be needed for hydrologic modeling and simulation when the preservation of both the seasonal and the annual statistics is desired (Salas and

The author used Univariate Seasonal PARMA(1,1) model to generate the sufficient training dataset, and obtained two generated samples. They included the input nodes/variables including mean temperature (Tmean), maximum temperature (Tmax), minimum temperature (Tmin), mean dew point temperature (DPmean), minimum relative humidity (RHmin), mean relative humidity (RHmean), mean wind speed (WSmean), maximum wind speed (WSmax), and sunshine duration (SD) in mean values and the output nodes/variables including PE and

**3.2 Construction of Univariate Seasonal PARMA(p,q) model** 

**3.1 Univariate seasonal periodic autoregressive moving average model** 

**3. Stochastic model** 

Abdelmohsen, 1993).

FAO-56 PM ETo in total values, respectively. Furthermore, they were generated for 100 years (Short-term), 500 years (Mid-term), and 1000 years (Long-term), respectively. The author selected the second generated sample, and the first 50 years of the 100, 500, and 1000 years was abandoned to eliminate the biases, respectively. The parameters of Univariate Seasonal PARMA(1,1) model were determined using the method of approximate least square, respectively.

### **4. Support Vector Machine Neural Networks Model (SVM-NNM)**

SVM-NNM has found wide application in several areas including pattern recognition, regression, multimedia, bio-informatics and artificial intelligence. Very recently, SVM-NNM is gaining recognition in hydrology (Dibike et al., 2001; Khadam & Kaluarachchi, 2004). SVM-NNM implements the structural risk minimization principle which attempts to minimize an upper bound on the generalization error by striking a right balance between the training performance error and the capacity of machine. The solution of traditional neural networks models such as MLP-NNM may tend to fall into a local optimal solution, whereas global optimum solution is guaranteed for SVM-NNM (Haykin, 2009). SVM-NNM is a new kind of classifier that is motivated by two concepts. First, transforming data into a high-dimensional space can transform complex problems into simpler problems that can use linear discriminant functions. Second, SVM-NNM is motivated by the concept of training and using only those inputs that are near the decision surface since they provide the most information about the classification. The first step in SVM-NNM is transforming the data into a high-dimensional space. This is done using radial basis function (RBF) that places a Gaussian at each sample data. Thus, the feature space becomes as large as the number of sample data. RBF uses backpropagation to train a linear combination of the gaussians to produce the final result. SVM-NNM, however, uses the idea of large margin classifiers for the training performance. This decouples the capacity of the classifier from the input space and at the same time provides good generalization. This is an ideal combination for classification (Vapnik, 1992, 2000; Principe et al., 2000; Tripathi et al., 2006).

In this study, the basic ideas of SVM-NNM are reviewed. Consider the finite training pattern i i x ,y . where <sup>n</sup> <sup>i</sup> x = a sample value of the input vector x considering of N training patterns; and <sup>n</sup> yi = the corresponding value of the desired model output. A nonlinear transformation function ( ) is defined to map the input space to a higher dimension feature space, nh . According to Cover's theorem (Cover, 1965), a linear function, f( ) , could be formulated in the high dimensional feature space to look for a nonlinear relationship between inputs and outputs in the original input space. It can be written as the following equation (3).

$$\overline{\mathbf{y}} = \mathbf{f}(\mathbf{x}) = \mathbf{w}^T \boldsymbol{\phi}(\mathbf{x}) + \mathbf{b} \tag{3}$$

where y = the actual model output. The coefficient w and b are adjustable model parameters. In SVM-NNM, we aim at minimizing the empirical risk. It can be written as the following equation (4).

$$\mathbf{R}\_{\rm emp} = \frac{1}{\mathbf{N}} \sum\_{i=1}^{N} \left| \mathbf{y}\_i - \overline{\mathbf{y}}\_i \right|\_{\rm c} \tag{4}$$

where Remp = the empirical risk; and i <sup>i</sup> y y = the Vapnik's ε-insensitive loss function. Following regularization theory (Haykin, 2009), the parameters w and b are calculated by minimizing the cost function. It can be written as the following equation (5).

$$\log \mu\_{\varepsilon}(\mathbf{w}, \boldsymbol{\xi}, \boldsymbol{\xi}^{\star}) = \frac{1}{2} \mathbf{w}^{\top} \mathbf{w} + \mathbf{C} \sum\_{i=1}^{N} (\boldsymbol{\xi}\_{i} + \boldsymbol{\xi}\_{i}^{\star}) \tag{5}$$

subject to the constraints: 1) i i <sup>i</sup> y y ε ξ i = 1, 2, ... , N, 2) i i <sup>i</sup> y y ε ξ i = 1, 2, ... , N, and 3) i i <sup>ξ</sup> ,<sup>ξ</sup> <sup>0</sup> i = 1, 2, ... , N. where ψε (w,ξ,<sup>ξ</sup> ) = the cost function; i i <sup>ξ</sup> ,ξ = positive slack variables; and C = the cost constant. The first term of the cost function, which represents weight decay, is used to regularize weight sizes and to penalize large weights. This helps in improving generalization performance (Hush and Horne, 1993). The second term of the cost function, which represents penalty function, penalizes deviations of y from y larger than using Vapnik's ε-insensitive loss function. The cost constant C determines the amount up to which deviations from ε are tolerated. Deviations above ε are denoted by i ξ , whereas deviations below – ε are denoted by i <sup>ξ</sup> . The constrained quadratic optimization problem can be solved using the method of Lagrangian multipliers (Haykin, 2009). From this solution, the coefficient w can be written as the following equation (6).

$$\mathbf{w} = \sum\_{i=1}^{N} (\mathbf{a}\_i - \mathbf{a}\_i^\*) \phi(\mathbf{x}\_i) \tag{6}$$

where αi i ,α = the Lagrange multipliers, which are positive real constants. The data points corresponding to non-zero values for i i (α α ) are called support vectors. In SVM-NNM to calculate PE and FAO-56 PM ETo, there are several possibilities for the choice of kernel function, including linear, polynomial, sigmoid, splines and RBF. In this study, RBF is used to map the input data into higher dimensional feature space. RBF can be written as the following equation (7).

$$\mathbf{k}(\mathbf{x}, \mathbf{x}\_{\mathbf{j}}) = \Phi\_1 = \exp(-\mathbf{B}\_1 \mathbf{R}\_{\mathbf{j}}^2) = \exp\left(-\frac{\left\|\mathbf{x}\_{\mathbf{i}} - \mathbf{x}\_{\mathbf{j}}\right\|^2}{2\sigma^2}\right) \tag{7}$$

where i, j = the input layer and the hidden layer; K(x,x )j <sup>1</sup> = the inner product kernel function; 1 <sup>2</sup> <sup>1</sup> <sup>B</sup> <sup>2</sup><sup>σ</sup> , and has a constant value; and σ = the width/spread of RBF, which can be adjusted to control the expressivity of RBF. The function for the single node of the output layer which receives the calculated results of RBF can be written as the following equation (8).

$$\mathbf{G}\_{\mathbf{k}} = [\sum\_{j=1}^{N} (\mathbf{a}\_{j} - \mathbf{a}\_{j}^{\*}) \cdot \mathbf{K}(\mathbf{x}, \mathbf{x}\_{j})] + \mathbf{B} \tag{8}$$

where k = the output layer; Gk = the calculated value of the single output node; and B = the bias in the output layer. Equation (8), finally, takes the form of equation (9) and (10), which represents SVM-NNM for the modeling of PE and FAO-56 PM ETo.

where Remp = the empirical risk; and i <sup>i</sup> y y = the Vapnik's ε-insensitive loss function. Following regularization theory (Haykin, 2009), the parameters w and b are calculated by

> T ε i i

subject to the constraints: 1) i i <sup>i</sup> y y ε ξ i = 1, 2, ... , N, 2) i i <sup>i</sup> y y ε ξ i = 1, 2, ... , N,

variables; and C = the cost constant. The first term of the cost function, which represents weight decay, is used to regularize weight sizes and to penalize large weights. This helps in improving generalization performance (Hush and Horne, 1993). The second term of the cost function, which represents penalty function, penalizes deviations of y from y larger than

 using Vapnik's ε-insensitive loss function. The cost constant C determines the amount up to which deviations from ε are tolerated. Deviations above ε are denoted by i ξ , whereas deviations below – ε are denoted by i <sup>ξ</sup> . The constrained quadratic optimization problem can be solved using the method of Lagrangian multipliers (Haykin, 2009). From this

ii i

i j 2

<sup>2</sup><sup>σ</sup> , and has a constant value; and σ = the width/spread of RBF, which can be

w (α α ) (x ) 

where αi i ,α = the Lagrange multipliers, which are positive real constants. The data points

calculate PE and FAO-56 PM ETo, there are several possibilities for the choice of kernel function, including linear, polynomial, sigmoid, splines and RBF. In this study, RBF is used to map the input data into higher dimensional feature space. RBF can be written as the

j 1 1j 2

where i, j = the input layer and the hidden layer; K(x,x )j <sup>1</sup> = the inner product kernel

adjusted to control the expressivity of RBF. The function for the single node of the output layer which receives the calculated results of RBF can be written as the following equation (8).

k jj j

G [(α α ) K(x,x )] B 

where k = the output layer; Gk = the calculated value of the single output node; and B = the bias in the output layer. Equation (8), finally, takes the form of equation (9) and (10), which

N

j 1

represents SVM-NNM for the modeling of PE and FAO-56 PM ETo.

k(x,x ) exp(-B R ) exp <sup>2</sup><sup>σ</sup>

1 <sup>ψ</sup> (w,ξ,<sup>ξ</sup> ) ww C (ξ ξ ) <sup>2</sup> 

N

(5)

= the cost function; i i <sup>ξ</sup> ,ξ = positive slack

(6)

(7)

are called support vectors. In SVM-NNM to

2

x x

(8)

i 1

minimizing the cost function. It can be written as the following equation (5).

solution, the coefficient w can be written as the following equation (6).

N

i 1

 

and 3) i i <sup>ξ</sup> ,<sup>ξ</sup> <sup>0</sup> i = 1, 2, ... , N. where ψε (w,ξ,<sup>ξ</sup> )

corresponding to non-zero values for i i (α α )

following equation (7).

function; 1 <sup>2</sup> <sup>1</sup> <sup>B</sup>

$$\text{PE} = \text{W}\_1 \cdot \Phi\_2(\text{G}\_{\text{k}}) = \text{W}\_1 \cdot \Phi\_2[[\sum\_{j=1}^{N} (\text{a}\_j - \text{a}\_j^\star) \cdot \text{K}(\text{x}, \text{x}\_j)] + \text{B}] \tag{9}$$

$$\text{FAO-56 PM ET}\_{\text{o}} = \text{W}\_{2} \cdot \Phi\_{2}(\text{G}\_{\text{k}}) = \text{W}\_{2} \cdot \Phi\_{2}[[\sum\_{j=1}^{N} (\text{a}\_{j} - \text{a}\_{j}^{\*}) \cdot \text{K}(\text{x}, \text{x}\_{j})] + \text{B}] \tag{10}$$

where Φ<sup>2</sup> ( ) = the linear sigmoid transfer function; W1 = the specific weights connected to the output variable of PE; and W2 = the specific weights connected to the output variable of FAO-56 PM ETo. A number of SVM-NNM computer programs are now available for the modeling of PE and FAO-56 PM ETo. NeuroSolutions 5.0 computer program was used to develop SVM-NNM structure. Fig. 2 shows the developed structure of SVM-NNM. From the Fig. 2, the input nodes/variables of climatic data are mean temperature (Tmean), maximum temperature (Tmax), minimum temperature (Tmin), mean dew point temperature (DPmean), minimum relative humidity (RHmin), mean relative humidity (RHmean), mean wind speed (WSmean), maximum wind speed (WSmax), and sunshine duration (SD) in mean values (01/1985-12/1990). The output nodes/variables of climatic data are PE and FAO-56 PM ETo in total values (01/1985-12/1990).

Fig. 2. The developed structure of SVM-NNM

#### **5. Study scope and data**

In this study, Gwangju and Haenam stations from the Yeongsan River catchment are selected among the 71 weather stations including Jeju-do under the control of the Korea meteorological administration (KMA). They have possessed long-term climatic data dating back over at least 30 years. The Yeongsan River catchment covers an area of 3455 km2, and lies between latitudes 34.4°N and 35.2°N, and between longitudes 126.2°E and 127.0°E. Fig. 3 shows the Yeongsan River catchment including Gwangju and Haenam stations. The climatic data, which was necessary for the modeling of PE and FAO-56 PM ETo, were collected from the Internet homepage of water management information system (www.wamis.go.kr) and the Korea meteorological administration (www.kma.go.kr).

Fig. 3. The Yeongsan River Catchment including Gwangju and Haenam stations

## **6. SVM-NNM performance**

## **6.1 Performance statistics**

The performance of SVM-NNM to account for calculating the monthly PE and FAO-56 PM ETo was evaluated using a wide variety of standard statistics index. A total of 3 different standard statistics indexes were employed; the coefficient of correlation (CC), root mean square error (RMSE), and Nash-Sutcliffe coefficient (R2) (Nash & Sutcliffe, 1970; ASCE, 1993). Table 1 shows summary of the statistics index in this study. where i y (x) = the calculated PE and FAO-56 PM ETo (mm/month); yi(x) = the observed PE and FAO-56 PM ETo (mm/month); uy = mean of the calculated PE and FAO-56 PM ETo (mm/month); uy = mean of the observed PE and FAO-56 PM ETo (mm/month); and n = total number of the monthly PE and FAO-56 PM ETo considered. A model which is effective in the modeling of PE and FAO-56 PM ETo accurately, and efficient in capturing the complex relationship among the various inputs and output variables involved in a particular problem, is considered the best. CC, RMSE, and R2 statistics quantify the efficiency of SVM-NNM in capturing the extremely complex, dynamic, and nonlinear relationships (Kim, 2011).


Table 1. Summary of statistics indexes

#### **6.2 Data normalization**

358 Evapotranspiration – Remote Sensing and Modeling

back over at least 30 years. The Yeongsan River catchment covers an area of 3455 km2, and lies between latitudes 34.4°N and 35.2°N, and between longitudes 126.2°E and 127.0°E. Fig. 3 shows the Yeongsan River catchment including Gwangju and Haenam stations. The climatic data, which was necessary for the modeling of PE and FAO-56 PM ETo, were collected from the Internet homepage of water management information system

(www.wamis.go.kr) and the Korea meteorological administration (www.kma.go.kr).

Fig. 3. The Yeongsan River Catchment including Gwangju and Haenam stations

The performance of SVM-NNM to account for calculating the monthly PE and FAO-56 PM ETo was evaluated using a wide variety of standard statistics index. A total of 3 different standard statistics indexes were employed; the coefficient of correlation (CC), root mean square error (RMSE), and Nash-Sutcliffe coefficient (R2) (Nash & Sutcliffe, 1970; ASCE, 1993). Table 1 shows summary of the statistics index in this study. where i y (x) = the calculated PE and FAO-56 PM ETo (mm/month); yi(x) = the observed PE and FAO-56 PM ETo (mm/month); uy = mean of the calculated PE and FAO-56 PM ETo (mm/month); uy = mean of the observed PE and FAO-56 PM ETo (mm/month); and n = total number of the monthly PE and FAO-56 PM ETo considered. A model which is effective in the modeling of PE and FAO-56 PM ETo accurately, and efficient in capturing the complex relationship among the various inputs and output variables involved in a particular problem, is considered the best. CC, RMSE, and R2 statistics quantify the efficiency of SVM-NNM in

capturing the extremely complex, dynamic, and nonlinear relationships (Kim, 2011).

**6. SVM-NNM performance 6.1 Performance statistics** 

The climatic data used in this study including mean temperature (Tmean), maximum temperature (Tmax), minimum temperature (Tmin), mean dew point temperature (DPmean), minimum relative humidity (RHmin), mean relative humidity (RHmean), mean wind speed (WSmean), maximum wind speed (WSmax), and sunshine duration (SD) were normalized for preventing and overcoming problem associated with the extreme values. An important reason for the normalization of input nodes is that each of input nodes represents an observed value in a different unit. Such input nodes are normalized, and the input nodes in dimensionless unit are relocated. The similarity effect of input nodes is thus eliminated (Kim et al., 2009). According to Zanetti et al. (2007), by grouping the daily values into averages, ETo may be calculated due to their highest stabilization. For data normalization, the data of input and output nodes were scaled in the range of [0 1] using the following equation (11).

$$\mathbf{Y}\_{\text{norm}} = \frac{\mathbf{Y}\_{\text{i}} - \mathbf{Y}\_{\text{min}}}{\mathbf{Y}\_{\text{max}} - \mathbf{Y}\_{\text{min}}} \tag{11}$$

where Ynorm = the normalized dimensionless data of the specific input node/variable; Yi = the observed data of the specific input node/variable; Ymin = the minimum data of the specific input node/variable; and Ymax = the maximum data of the specific input node/variable.

#### **6.3 Training performance**

The method for calculating parameters is generally called the training performance in the neural networks model category. The training performance of neural networks model is iterated until the training error is reached to the training tolerance. Iteration means one completely pass through a set of inputs and target patterns or data. In general, it is assumed that the neural networks model does not have any prior knowledge about the example problem before it is trained (Kim, 2004). A difficult task with the neural networks model is to choose the number of hidden nodes. The network geometry is problem dependent. This study adopted one hidden layer for the construction of SVM-NNM since it is well known that one hidden layer is enough to represent PE and FAO-56 PM ETo nonlinear complex relationship (Kumar et al., 2002; Zanetti et al., 2007). The testing performance in the modeling of PE and FAO-56 PM ETo, therefore, is carried out using the optimal parameters, which are calculated during the training performance.

The hybrid method, which was developed in this study, consisted of the following training patterns. First, the stochastic model was selected. As explained previously, Univariate Seasonal PARMA(1,1) model, which consisted of 1 pattern only, was used to generate the training dataset. Second, the data, which were generated by Univariate Seasonal PARMA(1,1) model, consisted of 3 patterns including 100 years (Short-term), 500 years (Mid-term), and 1000 years (Long-term), respectively. Finally, the neural networks model, which consisted of 1 pattern only including SVM-NNM, was used for the training and testing performances, respectively. Therefore, the hybrid method consisted of 3 training patterns including 100/PARMA(1,1)/SVM-NNM, 500/PARMA(1,1)/SVM-NNM, and 1000 /PARMA(1,1)/SVM-NNM, respectively. For Gwangju and Haenam stations, the training dataset including the climatic, PE, and FAO-56 PM ETo data were generated by Univariate Seasonal PARMA(1,1) model using observed data (01/1985-12/1990) for 100 years (Shortterm), 500 years (Mid-term), and 1000 years (Long-term), respectively. After the first 50 years of the generated data for 100, 500, and 1000 years was abandoned to eliminate the biases, the training performance should be carried out using SVM-NNM. Therefore, the total amount of data used for the training performance consisted of 600, 5400, and 11400, respectively. For the training performance of SVM-NNM, NeuroSolutions 5.0 computer program was used to carry out the training performance. Fig. 4(a)-(b) show SVM-NNM training performance using NeuroSolutions 5.0 computer program.


Fig. 4. SVM-NNM training performance using NeuroSolution 5.0 Computer Program

Table 2 shows the statistics results of the training performances for 3 training patterns of PE modeling. In PE of Gwangju station, from the Table 2, 100/PARMA(1,1)/SVM-NNM

that the neural networks model does not have any prior knowledge about the example problem before it is trained (Kim, 2004). A difficult task with the neural networks model is to choose the number of hidden nodes. The network geometry is problem dependent. This study adopted one hidden layer for the construction of SVM-NNM since it is well known that one hidden layer is enough to represent PE and FAO-56 PM ETo nonlinear complex relationship (Kumar et al., 2002; Zanetti et al., 2007). The testing performance in the modeling of PE and FAO-56 PM ETo, therefore, is carried out using the optimal parameters,

The hybrid method, which was developed in this study, consisted of the following training patterns. First, the stochastic model was selected. As explained previously, Univariate Seasonal PARMA(1,1) model, which consisted of 1 pattern only, was used to generate the training dataset. Second, the data, which were generated by Univariate Seasonal PARMA(1,1) model, consisted of 3 patterns including 100 years (Short-term), 500 years (Mid-term), and 1000 years (Long-term), respectively. Finally, the neural networks model, which consisted of 1 pattern only including SVM-NNM, was used for the training and testing performances, respectively. Therefore, the hybrid method consisted of 3 training patterns including 100/PARMA(1,1)/SVM-NNM, 500/PARMA(1,1)/SVM-NNM, and 1000 /PARMA(1,1)/SVM-NNM, respectively. For Gwangju and Haenam stations, the training dataset including the climatic, PE, and FAO-56 PM ETo data were generated by Univariate Seasonal PARMA(1,1) model using observed data (01/1985-12/1990) for 100 years (Shortterm), 500 years (Mid-term), and 1000 years (Long-term), respectively. After the first 50 years of the generated data for 100, 500, and 1000 years was abandoned to eliminate the biases, the training performance should be carried out using SVM-NNM. Therefore, the total amount of data used for the training performance consisted of 600, 5400, and 11400, respectively. For the training performance of SVM-NNM, NeuroSolutions 5.0 computer program was used to carry out the training performance. Fig. 4(a)-(b) show SVM-NNM

(a) Training Data (b) Training Performance

Fig. 4. SVM-NNM training performance using NeuroSolution 5.0 Computer Program

Table 2 shows the statistics results of the training performances for 3 training patterns of PE modeling. In PE of Gwangju station, from the Table 2, 100/PARMA(1,1)/SVM-NNM

which are calculated during the training performance.

training performance using NeuroSolutions 5.0 computer program.

training pattern produced the statistics results with CC value of 0.929, RMSE value of 16.360 mm and R2 value of 0.858. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.919, RMSE value of 17.942 mm and R2 value of 0.833. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.905, RMSE value of 20.489 mm and R2 value of 0.781, respectively. In PE of Haenam station, from the Table 2, 100/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.965, RMSE value of 9.643 mm and R2 value of 0.930. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.958, RMSE value of 10.673 mm and R2 value of 0.914. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.962, RMSE value of 10.105 mm and R2 value of 0.922, respectively. From the above results, the statistics results of the training performance for 100/PARMA(1,1)/SVM-NNM training pattern were better than those of the training performances for 500/PARMA(1,1)/SVM-NNM and 1000/ PARMA(1,1)/SVM-NNM training patterns for PE of Gwangju and Haenam stations, respectively.


Table 2. Statistics results of the training performances (PE)

Table 3 shows the statistics results of the training performances for 3 training patterns of FAO-56 PM ETo modeling. In FAO-56 PM ETo of Gwangju station, from the Table 3, 100/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.975, RMSE value of 8.517 mm and R2 value of 0.948. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.966, RMSE value of 10.237 mm and R2 value of 0.924. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.963, RMSE value of 11.061 mm and R2 value of 0.911, respectively. In FAO-56 PM ETo of Haenam station, from the Table 3, 100/PARMA(1,1) /SVM-NNM training pattern produced the statistics results with CC value of 0.965, RMSE value of 9.935 mm and R2 value of 0.926. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.956, RMSE value of 10.822 mm and R2 value of 0.912. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.962, RMSE value of 10.014 mm and R2 value of 0.925, respectively. From the above results, the statistics results of the training performance for 100/ PARMA(1,1)/SVM-NNM training pattern were better than those of the training performances for 500/PARMA(1,1)/SVM-NNM and 1000/PARMA(1,1)/SVM-NNM training patterns for FAO-56 PM ETo of Gwangju and Haenam stations, respectively. Therefore, the data length has less effect on the training performance of PE and FAO-56 PM ETo in this study. Tokar and Johnson (1999) suggested that the data length has less effect on the neural networks model performance than the data quality. Sivakumar et al. (2002) suggested that it is imperative to select a good training dataset from the available data series. They indicated that the best way to achieve a good training performance seems to be to include most of the extreme events such as very high and very low values in the training dataset. Furthermore, Kim (2011) did not carry out the statistics analysis of the training performance since the training dataset of 6 training patterns consisted of the generated (not observed) data only.


Table 3. Statistics results of the training performances (FAO-56 PM ETo)

## **6.4 Testing performance**

The neural networks model is tested by determining whether the model meets the objectives of modeling within some preestablished criteria or not. Of course, the optimal parameters, which are calculated during the training performance, are applied for the testing performance of the neural networks model (Kim, 2004). For the testing performance, the monthly climatic data (01/1985-12/1990) in Gwangju and Haenam stations were used. The total amount of data used for the testing performance consisted of 72 data for the monthly time series. The testing performance applied the cross-validation method in order to overcome the over-fitting problem of SVM-NNM. The cross-validation method is not to train all the training data until SVM-NNM reaches the minimum RMSE, but is to cross-validate with the testing data at the end of each training performance. If the over-fitting problem occurs, the convergence process over the mean square error of the testing data will not decrease but will increase as the training data are still trained (Bishop, 1994; Haykin, 2009).

Table 4 shows the statistics results of the testing performances for 3 training patterns of PE modeling. In PE of Gwangju station, from the Table 4, 100/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.955, RMSE value of 12.239 mm and R2 value of 0.908. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.956, RMSE value of 13.501 mm and R2 value of 0.888. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.953, RMSE value of 15.103 mm and R2 value of 0.860, respectively. In PE of Haenam station, from the Table 4, 100/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.966, RMSE value of 9.581 mm and R2 value of 0.932. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.968, RMSE value of 9.370 mm and R2 value of 0.935. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.953, RMSE value of 11.313 mm and R2 value of 0.905, respectively. From the above results, the statistics results of the

the neural networks model performance than the data quality. Sivakumar et al. (2002) suggested that it is imperative to select a good training dataset from the available data series. They indicated that the best way to achieve a good training performance seems to be to include most of the extreme events such as very high and very low values in the training dataset. Furthermore, Kim (2011) did not carry out the statistics analysis of the training performance since the training dataset of 6 training patterns consisted of the generated (not

> 500/PARMA(1,1)/ SVM-NNM

CC 0.975 0.966 0.963 RMSE (mm) 8.517 10.237 11.061 R2 0.948 0.924 0.911

CC 0.965 0.956 0.962 RMSE (mm) 9.935 10.822 10.014 R2 0.926 0.912 0.925

The neural networks model is tested by determining whether the model meets the objectives of modeling within some preestablished criteria or not. Of course, the optimal parameters, which are calculated during the training performance, are applied for the testing performance of the neural networks model (Kim, 2004). For the testing performance, the monthly climatic data (01/1985-12/1990) in Gwangju and Haenam stations were used. The total amount of data used for the testing performance consisted of 72 data for the monthly time series. The testing performance applied the cross-validation method in order to overcome the over-fitting problem of SVM-NNM. The cross-validation method is not to train all the training data until SVM-NNM reaches the minimum RMSE, but is to cross-validate with the testing data at the end of each training performance. If the over-fitting problem occurs, the convergence process over the mean square error of the testing data will not decrease but will increase as the training

Table 4 shows the statistics results of the testing performances for 3 training patterns of PE modeling. In PE of Gwangju station, from the Table 4, 100/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.955, RMSE value of 12.239 mm and R2 value of 0.908. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.956, RMSE value of 13.501 mm and R2 value of 0.888. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.953, RMSE value of 15.103 mm and R2 value of 0.860, respectively. In PE of Haenam station, from the Table 4, 100/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.966, RMSE value of 9.581 mm and R2 value of 0.932. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.968, RMSE value of 9.370 mm and R2 value of 0.935. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.953, RMSE value of 11.313 mm and R2 value of 0.905, respectively. From the above results, the statistics results of the

1000/PARMA(1,1) SVM-NNM

100/PARMA(1,1)/ SVM-NNM

Table 3. Statistics results of the training performances (FAO-56 PM ETo)

observed) data only.

Gwangju

 Haenam

Station Statistics

**6.4 Testing performance** 

Indexes

data are still trained (Bishop, 1994; Haykin, 2009).

testing performance for 100/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 500/PARMA(1,1)/SVM-NNM and 1000/ PARMA(1,1)/SVM-NNM training patterns for PE of Gwangju station. And, the statistics results of the testing performance for 500/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 100/PARMA(1,1)/SVM-NNM and 1000/ PARMA(1,1)/SVM-NNM training patterns for PE of Haenam station, respectively.



Table 5 shows the statistics results of the testing performances for 3 training patterns of FAO-56 PM ETo modeling. In FAO-56 PM ETo of Gwangju station, from the Table 5, 100/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.981, RMSE value of 7.300 mm and R2 value of 0.962. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.974, RMSE value of 8.803 mm and R2 value of 0.944. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.976, RMSE value of 9.448 mm and R2 value of 0.936, respectively. In FAO-56 PM ETo of Haenam station, from the Table 5, 100/PARMA(1,1)/ SVM-NNM training pattern produced the statistics results with CC value of 0.971, RMSE value of 9.007 mm and R2 value of 0.939. 500/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.970, RMSE value of 8.882 mm and R2 value of 0.941. 1000/PARMA(1,1)/SVM-NNM training pattern produced the statistics results with CC value of 0.972, RMSE value of 8.581 mm and R2 value of 0.945, respectively. From the above results, the statistics results of the testing performance for 100/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 500/PARMA(1,1)/SVM-NNM and 1000/PARMA(1,1)/SVM-NNM training patterns for FAO-56 PM ETo of Gwangju station. And, the statistics results of the testing performance for 1000/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 100/PARMA(1,1)/SVM-NNM and 500/PARMA(1,1)/SVM-NNM training patterns for FAO-56 PM ETo of Haenam station, respectively. Kim (2011) suggested, however, that the statistics results of testing performance for 1000/ PARMA(1,1)/GRNNM-GA training pattern were better than those of testing performances for 100/PARMA(1,1)/GRNNM-GA and 500/PARMA(1,1)/GRNNM-GA training patterns for the modeling of PE and ETr, South Korea. The continuous research will be needed to establish the neural networks models available for the optimal training patterns and modeling of PE and FAO-56 PM ETo.

From the statistics results of the testing performances for PE and FAO-56 PM ETo, the statistics results of FAO-56 PM ETo were better than those of PE. PE is the observed data as a reason and represents the natural phenomenon including strong nonlinear patterns and various uncertainties, whereas ETo is calculated by FAO-56 PM equation with the constant operation processes. In FAO-56 PM ETo, furthermore, the strong nonlinear patterns of the natural phenomena are transformed into linear patterns including the constant uncertainty. The author can consider that the modeling of FAO-56 PM ETo has significantly less uncertainty compared to that of PE. Kim (2011) suggested the similar results for the modeling of PE and ETr using the neural networks models. He suggested that the statistics results of ETr were better than those of PE for the modeling of PE and ETr using GRNNM-BP and GRNNM-GA, South Korea. Fig. 5(a)-(f) show comparison plots of observed and calculated PE/FAO-56 PM ETo for the testing performances of 3 training patterns for Gwangju station, respectively. Fig. 6(a)-(f) show comparison plots of observed and calculated PE/FAO-56 PM ETo for the testing performances of 3 training patterns for Haenam station, respectively.


Table 5. Statistics results of the testing performances (FAO-56 PM ETo)

### **6.5 Homogeneity evaluation**

FAO-56 PM ETo equation, which has been unanimously accepted by the FAO consultation members for ETo calculation (Allen et al., 1998), was used to calculate ETo since there are no observed data for ETo using a lysimeter, South Korea. Even if FAO-56 PM ETo is not observed data, the reliability for the calculated FAO-56 PM ETo is adequate and proper. Homogeneity evaluation was performed to compare the observed PE/FAO-56 PM ETo with the calculated PE/FAO-56 PM ETo for the results of the test performances, respectively. In this study, homogeneity evaluation consisted of the One-way ANOVA and the Mann-Whitney *U* test, respectively.

## **6.5.1 The One-way ANOVA**

The One-way ANOVA is a class of statistical analysis that is widely used because it encourages systematic decision making for the underlying problems that involve considerable uncertainty. It enables inferences to be made in such a way that sample data can be combined with statistical theory. It supposedly removes the effects of the biases of the individual, which leads to more rational and accurate decision making. The One-way ANOVA is the formal procedure for using statistical concepts and measures in performing decision making. The following six steps can be used to make statistical analysis of the Oneway ANOVA on the means and variances: 1) Formulation of hypotheses 2) Define the test statistic and its distribution 3) Specify the level of significance 4) Collect data and compute

reason and represents the natural phenomenon including strong nonlinear patterns and various uncertainties, whereas ETo is calculated by FAO-56 PM equation with the constant operation processes. In FAO-56 PM ETo, furthermore, the strong nonlinear patterns of the natural phenomena are transformed into linear patterns including the constant uncertainty. The author can consider that the modeling of FAO-56 PM ETo has significantly less uncertainty compared to that of PE. Kim (2011) suggested the similar results for the modeling of PE and ETr using the neural networks models. He suggested that the statistics results of ETr were better than those of PE for the modeling of PE and ETr using GRNNM-BP and GRNNM-GA, South Korea. Fig. 5(a)-(f) show comparison plots of observed and calculated PE/FAO-56 PM ETo for the testing performances of 3 training patterns for Gwangju station, respectively. Fig. 6(a)-(f) show comparison plots of observed and calculated PE/FAO-56 PM ETo for the testing performances of 3 training patterns for

> 100/PARMA(1,1)/ SVM-NNM

Table 5. Statistics results of the testing performances (FAO-56 PM ETo)

500/PARMA(1,1)/ SVM-NNM

CC 0.981 0.974 0.976 RMSE (mm) 7.300 8.803 9.448 R2 0.962 0.944 0.936

CC 0.971 0.970 0.972 RMSE (mm) 9.007 8.882 8.581 R2 0.939 0.941 0.945

FAO-56 PM ETo equation, which has been unanimously accepted by the FAO consultation members for ETo calculation (Allen et al., 1998), was used to calculate ETo since there are no observed data for ETo using a lysimeter, South Korea. Even if FAO-56 PM ETo is not observed data, the reliability for the calculated FAO-56 PM ETo is adequate and proper. Homogeneity evaluation was performed to compare the observed PE/FAO-56 PM ETo with the calculated PE/FAO-56 PM ETo for the results of the test performances, respectively. In this study, homogeneity evaluation consisted of the One-way ANOVA and the Mann-

The One-way ANOVA is a class of statistical analysis that is widely used because it encourages systematic decision making for the underlying problems that involve considerable uncertainty. It enables inferences to be made in such a way that sample data can be combined with statistical theory. It supposedly removes the effects of the biases of the individual, which leads to more rational and accurate decision making. The One-way ANOVA is the formal procedure for using statistical concepts and measures in performing decision making. The following six steps can be used to make statistical analysis of the Oneway ANOVA on the means and variances: 1) Formulation of hypotheses 2) Define the test statistic and its distribution 3) Specify the level of significance 4) Collect data and compute

1000/PARMA(1,1) SVM-NNM

Haenam station, respectively.

Station Statistics

**6.5 Homogeneity evaluation** 

Whitney *U* test, respectively.

**6.5.1 The One-way ANOVA** 

Gwangju

 Haenam Indexes

test statistic 5) Determine the critical value of the test statistic 6) Make a decision (McCuen, 1993; Salas et al., 2001; Ayyub and McCuen, 2003).

Fig. 5. Comparison plots of observed and calculated PE and FAO-56 PM ETo (Gwangju)

Fig. 6. Comparison plots of observed and calculated PE and FAO-56 PM ETo (Haenam)

(a) PE (100 year) (b) FAO-56 PM ETo (100 year)

(c) PE (500 year) (d) FAO-56 PM ETo (500 year)

(e) PE (1000 year) (f) FAO-56 PM ETo (1000 year)

Fig. 6. Comparison plots of observed and calculated PE and FAO-56 PM ETo (Haenam)

The One-way ANOVA on the means was performed and computed *t* statistic using twosample *t* test between the observed PE/FAO-56 PM ETo and the calculated PE/FAO-56 PM ETo, respectively. The critical value of *t* statistic was computed for the level of significance 5 percent (5%) and 1 percent (1%). If the computed value of *t* statistic is greater than the critical value of *t* statistic, the null hypothesis, which is the means are equal, should be rejected and the alternative hypothesis should be accepted. The One-way ANOVA on the variances was performed and computed *F* statistic using two-sample *F* test between the observed PE/FAO-56 PM ETo and the calculated PE/FAO-56 PM ETo, respectively. The critical value of *F* statistic was computed for the level of significance 5 percent (5%) and 1 percent (1%). If the computed value of *F* statistic is greater than the critical value of *F* statistic, the null hypothesis, which is the population variances are equal, should be rejected and the alternative hypothesis should be accepted.

Table 6 shows the results of the One-way ANOVA on the means of PE. The critical value of *t* statistic was computed as t0.05 =1.981 and t0.01=2.620 for the level of significance 5 percent (5%) and 1 percent (1%), respectively. The computed values of *t* statistic with 0.045 for 100/PARMA/SVM-NNM training pattern, 0.111 for 500/PARMA/SVM-NNM training pattern, 0.390 for 1000/PARMA/SVM-NNM training pattern were not significant for PE of Gwangju station. So, the null hypothesis, which is the means are equal, was accepted for PE of Gwangju station. Furthermore, the computed values of *t* statistic with 0.145 for 100/PARMA/SVM-NNM training pattern, 0.169 for 500/PARMA/SVM-NNM training pattern, 0.103 for 1000/PARMA/SVM-NNM training pattern were not significant for PE of Haenam station. So, the null hypothesis, which is the means are equal, was accepted for PE of Haenam station.


Table 6. Results of the One-way ANOVA on the means of PE

Table 7 shows the results of the One-way ANOVA on the means of FAO-56 PM ETo. The critical value of *t* statistic was computed as t0.05 =1.981 and t0.01=2.620 for the level of significance 5 percent (5%) and 1 percent (1%), respectively. The computed values of *t*  statistic with 0.040 for 100/PARMA/SVM-NNM training pattern, 0.283 for 500/PARMA /SVM-NNM training pattern, 0.483 for 1000/PARMA/SVM-NNM training pattern were not significant for FAO-56 PM ETo of Gwangju station. So, the null hypothesis, which is the means are equal, was accepted for FAO-56 PM ETo of Gwangju station. Furthermore, the computed values of *t* statistic with 0.231 for 100/PARMA/SVM-NNM training pattern, 0.071 for 500/PARMA/SVM-NNM training pattern, 0.018 for 1000/PARMA/SVM-NNM training pattern were not significant for FAO-56 PM ETo of Haenam station. So, the null hypothesis, which is the means are equal, was accepted for FAO-56 PM ETo of Haenam station.


Table 7. Results of the One-way ANOVA on the means of FAO-56 PM ETo

Table 8 shows the results of the One-way ANOVA on the variances of PE. The critical value of *F* statistic was computed as F0.05 =1.981 and F0.01=2.620 for the level of significance 5 percent (5%) and 1 percent (1%), respectively. The computed values of *F* statistic with 1.040 for 100/PARMA/SVM-NNM training pattern, 1.249 for 500/PARMA/SVM-NNM training pattern, 1.343 for 1000/PARMA/SVM-NNM training pattern were not significant for PE of Gwangju station. So, the null hypothesis, which is the variances are equal, was accepted for PE of Gwangju station. Furthermore, the computed values of *F* statistic with 1.033 for 100/PARMA/SVM-NNM training pattern, 1.030 for 500/PARMA/SVM-NNM training pattern, 1.036 for 1000/PARMA/SVM-NNM training pattern were not significant for PE of Haenam station. So, the null hypothesis, which is the variances are equal, was accepted for PE of Haenam station.


Table 8. Results of the One-way ANOVA on the variances of PE

Table 9 shows the results of the One-way ANOVA on the variances of FAO-56 PM ETo. The critical value of *F* statistic was computed as F0.05 =1.981 and F0.01=2.620 for the level of significance 5 percent (5%) and 1 percent (1%), respectively. The computed values of *F*  statistic with 1.055 for 100/PARMA/SVM-NNM training pattern, 1.045 for 500/PARMA/ SVM-NNM training pattern, 1.154 for 1000/PARMA/SVM-NNM training pattern were not significant for FAO-56 PM ETo of Gwangju station. So, the null hypothesis, which is the variances are equal, was accepted for FAO-56 PM ETo of Gwangju station. Furthermore, the computed values of *F* statistic with 1.033 for 100/PARMA/SVM-NNM training pattern, 1.021 for 500/PARMA/SVM-NNM training pattern, 1.031 for 1000/PARMA/SVM-NNM training pattern were not significant for FAO-56 PM ETo of Haenam station. So, the null hypothesis, which is the variances are equal, was accepted for FAO-56 PM ETo of Haenam station.


Table 9. Results of the One-way ANOVA on the variances of FAO-56 PM ETo

### **6.5.2 The Mann-Whitney** *U* **test**

368 Evapotranspiration – Remote Sensing and Modeling

Level of Significance

Table 8 shows the results of the One-way ANOVA on the variances of PE. The critical value of *F* statistic was computed as F0.05 =1.981 and F0.01=2.620 for the level of significance 5 percent (5%) and 1 percent (1%), respectively. The computed values of *F* statistic with 1.040 for 100/PARMA/SVM-NNM training pattern, 1.249 for 500/PARMA/SVM-NNM training pattern, 1.343 for 1000/PARMA/SVM-NNM training pattern were not significant for PE of Gwangju station. So, the null hypothesis, which is the variances are equal, was accepted for PE of Gwangju station. Furthermore, the computed values of *F* statistic with 1.033 for 100/PARMA/SVM-NNM training pattern, 1.030 for 500/PARMA/SVM-NNM training pattern, 1.036 for 1000/PARMA/SVM-NNM training pattern were not significant for PE of Haenam station. So, the null hypothesis, which is the variances are equal, was accepted for

> Level of Significance

Table 9 shows the results of the One-way ANOVA on the variances of FAO-56 PM ETo. The critical value of *F* statistic was computed as F0.05 =1.981 and F0.01=2.620 for the level of significance 5 percent (5%) and 1 percent (1%), respectively. The computed values of *F*  statistic with 1.055 for 100/PARMA/SVM-NNM training pattern, 1.045 for 500/PARMA/ SVM-NNM training pattern, 1.154 for 1000/PARMA/SVM-NNM training pattern were not significant for FAO-56 PM ETo of Gwangju station. So, the null hypothesis, which is the variances are equal, was accepted for FAO-56 PM ETo of Gwangju station. Furthermore, the computed values of *F* statistic with 1.033 for 100/PARMA/SVM-NNM training pattern, 1.021 for 500/PARMA/SVM-NNM training pattern, 1.031 for 1000/PARMA/SVM-NNM training pattern were not significant for FAO-56 PM ETo of Haenam station. So, the null hypothesis, which is the variances are equal, was accepted

100/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 1.040 Accept/Accept 500/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 1.249 Accept/Accept 1000/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 1.343 Accept/Accept

100/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 1.033 Accept/Accept 500/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 1.030 Accept/Accept 1000/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 1.036 Accept/Accept

Critical *F* Statistic

Table 7. Results of the One-way ANOVA on the means of FAO-56 PM ETo

100/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 0.040 Accept/Accept 500/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 0.283 Accept/Accept 1000/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 0.483 Accept/Accept

100/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 0.231 Accept/Accept 500/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 0.071 Accept/Accept 1000/PARMA(1,1)/SVM-NNM 0.05/0.01 1.981/2.620 0.018 Accept/Accept

Two-sample *t* test

Two-sample *F* test

Null Hypothesis

Computed *F* Statistic

Null Hypothesis

Computed *t* Statistic

Critical *t* Statistic

Station Training

Gwangju

Haenam

PE of Haenam station.

Gwangju

Haenam

Station Training

Pattern

for FAO-56 PM ETo of Haenam station.

Table 8. Results of the One-way ANOVA on the variances of PE

Pattern

The Mann-Whitney *U* test is a nonparametric alternative to the two-sample *t* test for two independent samples and can be used to test whether two independent samples have been taken from the same population. It is the most powerful alternative to the two-sample *t* test. Therefore, when the assumptions of the two-sample *t* test are violated or are difficult to evaluate such as with small samples, the Mann-Whitney *U* test should be applied. The Mann-Whitney *U* test is to be used in the case of two independent samples, and the Kruskal-Wallis test is an extension of the Mann-Whitney *U* test for the case of more than two independent samples (McCuen, 1993; Salas et al., 2001; Ayyub and McCuen, 2003).

The Mann-Whitney *U* test was performed and computed *z* statistic between the observed PE/FAO-56 PM ETo and the calculated PE/FAO-56 PM ETo, respectively. The critical value of *z* statistic was computed for the level of significance 5 percent (5%) and 1 percent (1%). If the computed value of *z* statistic is greater than the critical value of *z* statistic, the null hypothesis, which is the two independent samples are from the same population, should be rejected and the alternative hypothesis should be accepted in this study.

Table 10 shows the results of the Mann-Whitney *U* test of PE. The critical value of *z* statistic was computed as z0.05 =1.645 and z0.01=2.327 for the level of significance 5 percent (5%) and 1 percent (1%), respectively. The computed values of *z* statistic with -0.196 for 100/PARMA/ SVM-NNM training pattern, -0.136 for 500/PARMA/SVM-NNM training pattern, -0.288 for 1000/PARMA/SVM-NNM training pattern were not significant for PE of Gwangju station. So, the null hypothesis, which is the two independent samples are from the same population, was accepted for PE of Gwangju station. Furthermore, the computed values of *z* statistic with -0.172 for 100/PARMA/SVM-NNM training pattern, - 0.124 for 500/PARMA/SVM-NNM training pattern, -0.076 for 1000/PARMA/SVM-NNM training pattern were not significant for PE of Haenam station. So, the null hypothesis, which is the two independent samples are from the same population, was accepted for PE of Haenam station.

Table 11 shows the results of the Mann-Whitney *U* test of FAO-56 PM ETo. The critical value of *z* statistic was computed as z0.05 =1.645 and z0.01=2.327 for the level of significance 5 percent (5%) and 1 percent (1%), respectively. The computed values of *z* statistic with -0.056 for 100/PARMA/SVM-NNM training pattern, -0.515 for 500/PARMA/SVM-NNM training pattern, -0.711 for 1000/PARMA/SVM-NNM training pattern were not significant for FAO-56 PM ETo of Gwangju station. So, the null hypothesis, which is the two independent samples are from the same population, was accepted for FAO-56 PM ETo of Gwangju station. Furthermore, the computed values of *z* statistic with -0.380 for 100/PARMA/SVM-NNM training pattern, -0.212 for 500/PARMA/SVM-NNM training pattern, -0.176 for 1000/PARMA/SVM-NNM training pattern were not significant for FAO-56 PM ETo of Haenam station. So, the null hypothesis, which is the two independent samples are from the same population, was accepted for FAO-56 PM ETo of Haenam station.


Table 10. Results of the Mann-Whitney *U* test of PE


Table 11. Results of the Mann-Whitney *U* test of FAO-56 PM ETo

## **7. Construction of the Bivariate Linear Regression Analysis Model**

The bivariate linear regression analysis model (BLRAM) was adopted to calculate FAO-56 PM ETo simply using the observed PE and compare the observed PE and the calculated FAO-56 PM ETo. The BLRAM is a conventional and universal model, which can calculate FAO-56 PM ETo simply using the observed PE for Gwangju and Haenam stations, respectively. The BLRAM consists of two variables; PE as the independent variable Xt, and FAO-56 PM ETo as the dependent variable Yt. The mathematical expression can be written as the following equation (12) (McCuen, 1993; Salas et al., 2001; Ayyub and McCuen, 2003).

$$\mathbf{Y}\_{\mathbf{t}} = \mathbf{b}\_0 + \mathbf{b}\_1 \cdot \mathbf{X}\_{\mathbf{t}} \tag{12}$$

where b1 = the slope coefficient, which is also known as the regression coefficient because it is calculated by the result of a regression analysis. Using the BLRAM, the correlation relationship was investigated between the observed PE and FAO-56 PM ETo for 3 training patterns. A very good relationship was found with the BLRAM, which could calculate FAO-56 PM ETo in this study. The results of the BLRAM were the same for 3 training patterns. Therefore, it can be considered that the observed PE and FAO-56 PM ETo for 3 training patterns are homogeneous groups. It can be inferred from the homogeneity evaluation of the previous chapter.

station. Furthermore, the computed values of *z* statistic with -0.380 for 100/PARMA/SVM-NNM training pattern, -0.212 for 500/PARMA/SVM-NNM training pattern, -0.176 for 1000/PARMA/SVM-NNM training pattern were not significant for FAO-56 PM ETo of Haenam station. So, the null hypothesis, which is the two independent samples are from the

> Level of Significance

> Level of Significance

100/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.196 Accept/Accept 500/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.136 Accept/Accept 1000/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.288 Accept/Accept

100/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.172 Accept/Accept 500/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.124 Accept/Accept 1000/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.076 Accept/Accept

100/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.056 Accept/Accept 500/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.515 Accept/Accept 1000/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.711 Accept/Accept

100/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.380 Accept/Accept 500/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.212 Accept/Accept 1000/PARMA(1,1)/SVM-NNM 0.05/0.01 1.645/2.327 -0.176 Accept/Accept

Critical *z* Statistic

Critical *z* Statistic

Y b bX t 0 1t (12)

Mann-Whitney *U* test

Mann-Whitney *U* test

Computed *z* Statistic

Null Hypothesis

Null Hypothesis

Computed *z* Statistic

same population, was accepted for FAO-56 PM ETo of Haenam station.

Station Training

Station Training

Gwangju

Haenam

Gwangju

Haenam

the previous chapter.

Pattern

Table 10. Results of the Mann-Whitney *U* test of PE

Pattern

Table 11. Results of the Mann-Whitney *U* test of FAO-56 PM ETo

**7. Construction of the Bivariate Linear Regression Analysis Model** 

The bivariate linear regression analysis model (BLRAM) was adopted to calculate FAO-56 PM ETo simply using the observed PE and compare the observed PE and the calculated FAO-56 PM ETo. The BLRAM is a conventional and universal model, which can calculate FAO-56 PM ETo simply using the observed PE for Gwangju and Haenam stations, respectively. The BLRAM consists of two variables; PE as the independent variable Xt, and FAO-56 PM ETo as the dependent variable Yt. The mathematical expression can be written as the following equation (12) (McCuen, 1993; Salas et al., 2001; Ayyub and McCuen, 2003).

where b1 = the slope coefficient, which is also known as the regression coefficient because it is calculated by the result of a regression analysis. Using the BLRAM, the correlation relationship was investigated between the observed PE and FAO-56 PM ETo for 3 training patterns. A very good relationship was found with the BLRAM, which could calculate FAO-56 PM ETo in this study. The results of the BLRAM were the same for 3 training patterns. Therefore, it can be considered that the observed PE and FAO-56 PM ETo for 3 training patterns are homogeneous groups. It can be inferred from the homogeneity evaluation of Table 12 shows the BLRAM, goodness-of-fit test, and regression coefficient between the observed PE and FAO-56 PM ETo for 3 training patterns of Gwangju and Haenam stations, respectively. From the Table 12, for Gwangju station, the regression coefficient of the BLRAM indicates that FAO-56 PM ETo increases 0.9060 mm as each 1 mm increase in PE. R2 = 0.966 indicates that the total variance of FAO-56 PM ETo corresponds to 96.6%. The ratio of S(e)/S(y) = 0.187 suggests a very good level of accuracy. In addition, the standard error of estimate (SEE) decreases from 37.471 mm (S(y)) to 7.007 mm (S(e)). where S(y) = the standard deviation of FAO-56 PM ETo; and S(e) = the standard error of FAO-56 PM ETo using the equation (12). The overall deviations are nearly zero, and this tendency always occurs for the BLRAM. The standard error ratios of the regression coefficient (b1) and the intercept (b0) are 0.023 and 0.911, which indicates that the regression coefficient is relatively more accurate than the intercept. For Haenam station, the regression coefficient of the BLRAM indicates that FAO-56 PM ETo increases 0.9574 mm as each 1 mm increase in PE. R2 = 0.919 indicates that the total variance of FAO-56 PM ETo corresponds to 91.9%. The ratio of S(e)/S(y) = 0.287 suggests a very good level of accuracy. In addition, the standard error of estimate (SEE) decreases from 36.794 mm (S(y)) to 10.560 mm (S(e)). where S(y) = the standard deviation of FAO-56 PM ETo; and S(e) = the standard error of FAO-56 PM ETo using the equation (12). The overall deviations are nearly zero, and this tendency always occurs for the BLRAM. The standard error ratios of the regression coefficient (b1) and the intercept (b0) are 0.035 and 0.345, which indicates that the regression coefficient is relatively more accurate than the intercept. If a large and negative intercept exists, it can create some problems for forecasting or modeling (McCuen, 1993). Fig. 7(a)-(b) show comparison plots of the observed PE and FAO-56 PM ETo for 1000/PARMA(1,1)/SVM-NNM training pattern of Gwangju and Haenam station, respectively.


Table 12. Regression analysis between the observed PE and FAO-56 PM ETo

Fig. 7. Comparison plots of the observed PE and FAO-56 PM ETo

## **8. Conclusions**

The hybrid method was developed for the modeling of the monthly PE and FAO-56 PM ETo simultaneously. The author determined in advance 4 kinds of Univariate Seasonal PARMA(p,q) models including PARMA(1,1), PARMA(1,2), PARMA (2,1), and PARMA(2,2), which are the low-order models and contain the seasonal properties. As a result, the author selected Univariate Seasonal PARMA(1,1) model, which show the best statistical properties and is simple in parameter estimation. The data which were generated by Univariate Seasonal PARMA(1,1) model consisted of 100 years (Short-term), 500 years (Mid-term), and 1000 years (Long-term), respectively. The following conclusions can be drawn from this study.

[1] The statistics results of the testing performance for 100/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 500/PARMA(1,1) /SVM-NNM and 1000/PARMA(1,1)/SVM-NNM training patterns for PE of Gwangju station. And, the statistics results of the testing performance for 500/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 100 /PARMA(1,1)/SVM-NNM and 1000/ PARMA(1,1)/SVM-NNM training patterns for PE of Haenam station, respectively

[2] The statistics results of the testing performance for 100/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 500/PARMA(1,1) /SVM-NNM and 1000/PARMA(1,1)/SVM-NNM training patterns for FAO-56 PM ETo of Gwangju station. And, the statistics results of the testing performance for 1000 /PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 100/PARMA(1,1)/SVM-NNM and 500/PARMA(1,1)/SVM-NNM training patterns for FAO-56 PM ETo of Haenam station, respectively

[3] Homogeneity evaluation consisted of the One-way ANOVA and the Mann-Whitney *U* test. The null hypothesis, which is the means are equal, was accepted using the One-way ANOVA on the means for PE and FAO-56 PM ETo of Gwangju and Haenam stations, respectively. And, the null hypothesis, which is the variances are equal, was accepted using the One-way ANOVA on the variances for PE and FAO-56 PM ETo of Gwangju and Haenam stations, respectively. The null hypothesis, which is the two independent samples are from the same population, was accepted using the Mann-Whitney *U* test for PE and FAO-56 PM ETo of Gwangju and Haenam stations, respectively.

[4] The BLRAM was adopted to calculate FAO-56 PM ETo simply using the observed PE and compare the observed PE and the calculated FAO-56 PM ETo of Gwangju and Haenam stations, respectively. A very good relationship was found with the BLRAM, which could calculate FAO-56 PM ETo.

As PE and FAO-56 PM ETo are relatively important for the design of irrigation facilities and agricultural reservoirs, the spread of an automatic measuring system for PE and FAO-56 PM ETo is important and urgent to ensure the reliable and accurate data from the measurements of PE and FAO-56 PM ETo. Furthermore, the continuous research will be needed to establish the neural networks models available for the optimal training patterns and modeling of PE and FAO-56 PM ETo.

## **9. References**

Allen, R.G.; Jensen, M.E.; Wright, J.L. & Burman, R.D. (1989). Operational estimates of reference evapotranspiration. *Agronomy Journal*, Vol. 81, No. 4, pp. 650-662, ISSN: 0002-1962.

The hybrid method was developed for the modeling of the monthly PE and FAO-56 PM ETo simultaneously. The author determined in advance 4 kinds of Univariate Seasonal PARMA(p,q) models including PARMA(1,1), PARMA(1,2), PARMA (2,1), and PARMA(2,2), which are the low-order models and contain the seasonal properties. As a result, the author selected Univariate Seasonal PARMA(1,1) model, which show the best statistical properties and is simple in parameter estimation. The data which were generated by Univariate Seasonal PARMA(1,1) model consisted of 100 years (Short-term), 500 years (Mid-term), and 1000 years (Long-term), respectively. The following conclusions can be drawn from this

[1] The statistics results of the testing performance for 100/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 500/PARMA(1,1) /SVM-NNM and 1000/PARMA(1,1)/SVM-NNM training patterns for PE of Gwangju station. And, the statistics results of the testing performance for 500/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 100 /PARMA(1,1)/SVM-NNM and 1000/ PARMA(1,1)/SVM-NNM training patterns for PE of

[2] The statistics results of the testing performance for 100/PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 500/PARMA(1,1) /SVM-NNM and 1000/PARMA(1,1)/SVM-NNM training patterns for FAO-56 PM ETo of Gwangju station. And, the statistics results of the testing performance for 1000 /PARMA(1,1)/SVM-NNM training pattern were better than those of the testing performances for 100/PARMA(1,1)/SVM-NNM and 500/PARMA(1,1)/SVM-NNM training

[3] Homogeneity evaluation consisted of the One-way ANOVA and the Mann-Whitney *U* test. The null hypothesis, which is the means are equal, was accepted using the One-way ANOVA on the means for PE and FAO-56 PM ETo of Gwangju and Haenam stations, respectively. And, the null hypothesis, which is the variances are equal, was accepted using the One-way ANOVA on the variances for PE and FAO-56 PM ETo of Gwangju and Haenam stations, respectively. The null hypothesis, which is the two independent samples are from the same population, was accepted using the Mann-Whitney *U* test for PE and

[4] The BLRAM was adopted to calculate FAO-56 PM ETo simply using the observed PE and compare the observed PE and the calculated FAO-56 PM ETo of Gwangju and Haenam stations, respectively. A very good relationship was found with the BLRAM, which could

As PE and FAO-56 PM ETo are relatively important for the design of irrigation facilities and agricultural reservoirs, the spread of an automatic measuring system for PE and FAO-56 PM ETo is important and urgent to ensure the reliable and accurate data from the measurements of PE and FAO-56 PM ETo. Furthermore, the continuous research will be needed to establish the neural networks models available for the optimal training patterns and modeling of PE

Allen, R.G.; Jensen, M.E.; Wright, J.L. & Burman, R.D. (1989). Operational estimates of

reference evapotranspiration. *Agronomy Journal*, Vol. 81, No. 4, pp. 650-662, ISSN:

patterns for FAO-56 PM ETo of Haenam station, respectively

FAO-56 PM ETo of Gwangju and Haenam stations, respectively.

**8. Conclusions** 

Haenam station, respectively

calculate FAO-56 PM ETo.

and FAO-56 PM ETo.

0002-1962.

**9. References** 

study.


Khadam, I.M. & Kaluarachchi, J.J. (2004). Use of soft information to describe the relative

Kim, S. (2004). Neural Networks Model and Embedded Stochastic Processes for

Kim, S. (2011). Nonlinear hydrologic modeling using the stochastic and neural networks approach. *Disaster Advances,* Vol. 4, No. 1, pp. 53-63, ISSN: 0974-262X. Kim, S.; Kim, J.H. & Park, K.B. (2009). Neural networks models for the flood forecasting and

Kim, S. & Kim, H.S. (2008). Neural networks and genetic algorithm approach for nonlinear

Kisi, O. (2006). Generalized regression neural networks for evapotranspiration modeling. *Hydrological Sciences Journal*, Vol. 51, No. 6, pp. 1092-1105, ISSN: 0262-6667. Kisi, O. (2007). Evapotranspiration modeling from climatic data using a neural computing technique. *Hydrological Processes,* Vol. 21, pp. 1925-1934, ISSN: 0885-6087. Kisi, O. (2009). Modeling monthly evaporation using two different neural computing technique. *Irrigation Science,* Vol. 27, No. 5, pp. 417-430, ISSN: 0342-7188. Kisi, O. & Ozturk, O. (2007). Adaptive neurofuzzy computing technique for

Kumar, M.; Bandyopadhyay, A.; Raghuwanshi, N.S. & Singh, R. (2008). Comparative study

Kumar, M.; Raghuwanshi, N.S. & Singh, R. (2009). Development and validation of GANN

Kumar, M.; Raghuwanshi, N.S.; Singh, R.; Wallender, W.W. & Pruitt, W.O. (2002).

McCuen, R.H. (1993). *Microcomputer applications in statistical hydrology,* Prentice Hall, ISBN:

McKenzie, R.S. & Craig, A.R. (2001). Evaluation of river losses from the Orange River using hydraulic modeling. *Journal of Hydrology,* Vol. 241, pp. 62-69, ISSN: 0022-1694. Mishra, A.K.; Desai, V.R. & Singh, V.P. (2007). Drought forecasting using a hybrid stochastic

Monteith, J.L. (1965). The state and movement of water in living organism, *Proceedings of* 

*Irrigation Science,* Vol. 26, No. 6, pp. 531-545, ISSN: 0342-7188.

*Water Management,* Vol. 95, No. 5, pp. 553-565, ISSN: 0378-3774.

40, No. 11, W11505, ISSN: 0043-1397.

1, pp. 141-148, ISSN: 1226-7988.

pp. 51-63, ISSN: 0974-262X.

299-317, ISSN: 0022-1694.

Vol. 133, No. 4, pp. 368-379, ISSN: 0733-9437.

Vol. 14, No. 2, pp. 131-140, ISSN: 1084-0699.

0135852900, Englewood Cliffs, NJ, USA.

pp. 626-638, ISSN: 1084-0699.

NY, USA.

uncertainty of calibration data in hydrologic models. *Water Resources Research,* Vol.

Hydrological Analysis in South Korea. *KSCE Journal of Civil Engineering,* Vol. 8, No.

disaster prevention system in the small catchment. *Disaster Advances,* Vol. 2, No. 3,

evaporation and evapotranspiration modeling. *Journal of Hydrology,* Vol. 351, pp.

evapotranspiration estimation. *Journal of Irrigation and Drainage Engineering, ASCE,*

of conventional and artificial neural network-based ETo estimation models.

model for evapotranspiration estimation. *Journal of Hydrologic Engineering, ASCE*,

Estimating evapotranspiration using artificial neural network. *Journal of Irrigation and Drainage Engineering, ASCE,* Vol. 128, No. 4, pp. 224-233, ISSN: 0733-9437. Landeras, G.; Ortiz-Barredo, A. & Lopez, J.J. (2008). Comparison of artificial neural network

models and empirical and semi-empirical equations for daily reference evapotranspiration estimation in the Basque country (Northern Spain). *Agricultural* 

and neural network model. *Journal of Hydrologic Engineering, ASCE*, Vol. 12, No. 6,

*Evaporation and Environment*, pp. 205-234, Swansea, Cambridge University Press,


## **Modelling Evapotranspiration and the Surface Energy Budget in Alpine Catchments**

Giacomo Bertoldi1, Riccardo Rigon2 and Ulrike Tappeiner1,3 <sup>1</sup> *Institute for Alpine Environment, EURAC Research, Bolzano.* <sup>2</sup> *University of Trento.* <sup>3</sup> *Institute of Ecology, University of Innsbruck.* 1,2 *Italy* <sup>3</sup> *Austria.*

## **1. Introduction**

376 Evapotranspiration – Remote Sensing and Modeling

Trajkovic, S. (2005). Temperature-based approaches for estimating reference

Trajkovic, S; Todorovic, B. & Stankovic, M. (2003). Forecasting reference evapotranspiration

Tripathi, S.; Srinivas, V.V. & Nanjundish, R.S. (2006). Downscaling of precipitation for

Vallet-Coulomb, C.; Legesse, D.; Gasse, F.; Travi, Y. & Chernet, T. (2001). Lake evaporation

Vapnik, V.N. (1992). Principles of risk minimization for learning theory. In: *Advances in* 

Vapnik, V.N. (2010). *The nature of statistical learning theory 2nd Edition,* Springer-Verlag, ISBN:

Zanetti, S.S.; Sousa, E.F.; Oliveira, V.P.S.; Almeida, F.T. & Bernardo, S. (2007). Estimating

No. 4, pp. 316-323, ISSN: 0733-9437.

Vol. 330, pp. 621-640, ISSN: 0022-1694.

pp. 1-18, ISSN: 0022-1694.

0387987800, NY, USA.

ISSN: 0733-9437.

Vol. 129, No. 6, pp. 454-457, ISSN: 0733-9437.

831-838, Elsevier, ISBN: 1558602224, NY, USA.

evapotranspiration. *Journal of Irrigation and Drainage Engineering, ASCE,* Vol. 131,

by artificial neural networks. *Journal of Irrigation and Drainage Engineering, ASCE,*

climate change scenarios: A support vector machine approach. *Journal of Hydrology*,

estimates in tropical Africa (Lake Ziway, Ethiopia). *Journal of Hydrology*, Vol. 245,

*Neural Information Processing Systems Vol. 4*, Moody, Hanson & Lippmann, (Ed.), pp.

evapotranspiration using artificial neural network and minimum climatological data. *Journal of Irrigation and Drainage Engineering, ASCE,* Vol. 133, No. 2, pp. 83-89,

Accurate modelling of evapotranspiration (ET) is required to predict the effects of climate and land use changes on water resources, agriculture and ecosystems. Significant progress has been made in estimating ET at the global and regional scale. However, further efforts are needed to improve spatial accuracy and modeling capabilities in alpine regions (Brooks & Vivoni, 2008b). This chapter will point out the components of the energy budget needed to model ET, to discuss the fundamental equations and to provide an extended review of the parametrizations available in the hydrological and land surface models literature. The second part of the chapter will explore the complexity of the energy budget with special reference to mountain environments.

## **2. The energy budget components**

Evapotranspiration is controlled by the surface water and energy budget. In this section the single components of the energy budget will be discussed: radiation, soil heat flux, sensible and latent heat fluxes.

The surface energy budget inside a control volume can be written as:

$$
\Delta t \left( R\_{\text{fl}} - H - LE - G \right) = \Delta E \tag{1}
$$

where the energy fluxes concerning the soil-atmosphere interface in the time interval Δ*t* are:


The control volume is assumed with a thickness of some meters, so as to include the soil layer close to the surface and the first meters of the atmosphere, including the possible vegetation cover.

Besides the surface energy budget also the mass budget must be considered in order to quantify ET, namely the water conservation inside the control volume.

$$\frac{\Delta S}{\Delta t} = P - ET - R - R\_G - R\_S \tag{2}$$

where Δ*S* is the storing of the various supplies (underground and surface storage, soil moisture, vegetation interception, storing in channels); *P* is precipitation; *ET* = *λLE* is evapotranspiration; *λ* is the latent heat of vaporization; *R* is the surface runoff; *RG* is the runoff towards the deep water table; *RS* is the sub-surface runoff.

In the next sections it will be explained how the different components of the energy budget are usually described in hydrological and Land Surface Models (LSMs): radiation, sensible and latent heat flux, soil heat flux.

#### **2.1 Radiation**

The radiation is usually divided in short-wave components - indicated here as *SW* (including visible light, part of the ultraviolet radiation and the close infrared) and long-wave components - indicated here as *LW* (infrared radiation), with wavelength *λ* ranging respectively from 0.1 to 3 *μm* (98% of extraterrestrial radiation) and from 3 to 100 *μm* (2% of extraterrestrial radiation). In photosynthetic processes the short-wave radiation is further divided in photosynthetically active radiation (0.4 < *λ* < 0.7*μm*) and near infrared (Bonan, 1996).

Moreover, there is a further distinction between diffuse radiation *D* (coming by diffusion in the atmosphere from any direction) and direct radiation *P* (coming only from the sun), and also between radiation from the sky downwards ("*down*" ↓) and radiation from the soil upwards ("*up*" ↑).

The net radiation at the soil level can be factorized as follows:

$$R\_H = sw \cdot \left(R \downarrow\_{SW \, P} + R \downarrow\_{LW \, P} \right) + V \cdot \left(R \downarrow\_{SW \, D} + R \downarrow\_{LW \, D} \right) - \tag{3}$$

$$R \uparrow\_{SW} - R \uparrow\_{LW \, R} - R \uparrow\_{LW} + R \downarrow\_{SW \, O} + R \downarrow\_{LW \, O}$$

with:


*R* ↓*SW O* +*R* ↓*LW O* reflected radiation emitted by the surrounding surfaces;

The effects of topography on the diffused radiation can be expressed through the sky view factor *V*, indicating the sky fraction visible in one point.

In the presence of reliefs only part of the horizon is visible and this has consequences on the radiative exchanges.

*V* is defined as:

$$V = \omega / 2\pi \tag{4}$$

Fig. 1. Scheme of the solar radiation components.

2 Will-be-set-by-IN-TECH

Besides the surface energy budget also the mass budget must be considered in order to

where Δ*S* is the storing of the various supplies (underground and surface storage, soil moisture, vegetation interception, storing in channels); *P* is precipitation; *ET* = *λLE* is evapotranspiration; *λ* is the latent heat of vaporization; *R* is the surface runoff; *RG* is the

In the next sections it will be explained how the different components of the energy budget are usually described in hydrological and Land Surface Models (LSMs): radiation, sensible

The radiation is usually divided in short-wave components - indicated here as *SW* (including visible light, part of the ultraviolet radiation and the close infrared) and long-wave components - indicated here as *LW* (infrared radiation), with wavelength *λ* ranging respectively from 0.1 to 3 *μm* (98% of extraterrestrial radiation) and from 3 to 100 *μm* (2% of extraterrestrial radiation). In photosynthetic processes the short-wave radiation is further divided in photosynthetically active radiation (0.4 < *λ* < 0.7*μm*) and near infrared (Bonan,

Moreover, there is a further distinction between diffuse radiation *D* (coming by diffusion in the atmosphere from any direction) and direct radiation *P* (coming only from the sun), and also between radiation from the sky downwards ("*down*" ↓) and radiation from the soil upwards

*R* ↓*SW O* +*R* ↓*LW O* reflected radiation emitted by the surrounding surfaces;

The effects of topography on the diffused radiation can be expressed through the sky view

In the presence of reliefs only part of the horizon is visible and this has consequences on the

*V* = *ω*/2*π* (4)

*Rn* = *sw* · (*R* ↓*SW P* +*R* ↓*LW P*) + *V* · (*R* ↓*SW D* +*R* ↓*LW D*) − (3) *R* ↑*SW* −*R* ↑*LW R* −*R* ↑*LW* +*R* ↓*SW O* +*R* ↓*LW O*

<sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>P</sup>* <sup>−</sup> *ET* <sup>−</sup> *<sup>R</sup>* <sup>−</sup> *RG* <sup>−</sup> *RS* (2)

quantify ET, namely the water conservation inside the control volume.

Δ*S*

runoff towards the deep water table; *RS* is the sub-surface runoff.

The net radiation at the soil level can be factorized as follows:

*R* ↓*SW P* short-wave direct radiation; *R* ↓*LW P* long-wave direct radiation; *R* ↓*SW D* short-wave diffuse radiation; *R* ↓*LW D* long-wave diffuse radiation; *R* ↑*SW* short-wave reflected radiation; *R* ↑*LW R* long-wave reflected radiation; *R* ↑*LW* radiation emitted by surface;

factor *V*, indicating the sky fraction visible in one point.

and latent heat flux, soil heat flux.

**2.1 Radiation**

1996).

("*up*" ↑).

with:

radiative exchanges. *V* is defined as:

where *ω* is the solid angle seen from the point considered.

The presence of shadows due to surrounding mountains can be expressed through a factor *sw*, a function of topography and sun position, defined as:

$$sw = \begin{cases} 1 \text{ if the point is in the sun} \\ 0 \text{ if the point is in the shadow} \end{cases} \tag{5}$$

All direct radiation terms have to be multiplied by this factor.

In the next paragraphs we analyze in detail the parametrization of the single terms composing the radiation flux.

## **2.1.1 Direct radiation** *R* ↓*SW P* **and** *R* ↓*LW P*

The direct long-wave radiation *R* ↓*LW P* is emitted directly by the sun and therefore it is negligible at the soil level (differently from the long-wave diffuse radiation).

Usually the short-wave radiation *R* ↓*SW P* is assumed as an input variable, measured or calculated by an atmospheric model. The direct radiation can be written as the product of the extraterrestrial radiation *RExtr* by an attenuation factor varying in time and space.

$$R \downarrow\_{SW} P = F\_{att} R\_{Extr} \tag{6}$$

The extraterrestrial radiation can be easily calculated on the basis of geometric formulas (Iqbal, 1983). The atmospheric attenuation is due to Rayleigh diffusion, to the absorption on behalf of ozone and water vapor, to the extinction (both diffusion and absorption) due to atmospheric dust and shielding caused by the possible cloud cover. Moreover the absorption entity depends on the ray path length through the atmosphere, a function of the incidence angle and of the measurement point elevation. The effect of the latter can be very important in a mountain environment, where it is necessary to consider the shading effects.

Part of the dispersed radiation is then returned as short-wave diffuse radiation (*R* ↓*SW D*) and part of the energy absorbed by atmosphere is then re-emitted as long-wave diffuse radiation (*R* ↓*LW D*).

From a practical point of view, according to the application type and depending on the measured data possessed, the attenuation coefficient can be calculated with different degrees of complexity. The radiation transfer through the atmosphere is a well studied phenomenon and there exist many models providing the soil incident radiation spectrum in a detailed way, considering the various attenuation effects separately (Kondratyev, 1969).

## **2.1.2 Diffuse downward short-wave radiation** *R* ↓*SW D*

This term is a function of the atmospheric radiation due to Rayleigh dispersion and to the aereosols dispersion, as well as to the presence of cloud cover. The *R* ↓*SW D* actually is not isotropic and it depends on the sun position above the horizon. For its parametrization, see, for example, Paltrige & Platt (1976).

## **2.1.3 Diffuse downward long-wave radiation** *R* ↓*LW D*

Often this term in not provided by standard meteorological measurements, and many LSMs provide expressions to calculate it. This term indicates the long-wave radiation emitted by atmosphere towards the earth. It can be calculated starting from the knowledge of the distribution of temperature, humidity and carbon dioxide of the air column above. If this information is not available, various formulas, based only on ground measurements, can be found in literature with expressions as follows:

$$R \downarrow\_{LW} \mathbf{D} = \varepsilon\_a \sigma T\_a^4 \tag{7}$$

with:

*Ta* air temperature [K];

*�<sup>a</sup>* atmosphere emissivity *f*(*ea*, *Ta*, cloud cover);

*ea* vapor pressure [mb];

Usually for *�<sup>a</sup>* empirical formulas have been used, but it is also possible to provide a derivation based on physical topics like in Prata (1996). The cloud cover effect on this term is significant and not easy to consider in a simple way. Cloud cover data can be provided during the day by ground or satellite observations but, especially on night, is difficult to collect.

## **2.1.4 Reflected short-wave radiation** *R* ↑*SW*

This term indicates the short-wave energy reflection.

$$R \uparrow\_{SW} = a(R \downarrow\_{SW \mid P} + R \downarrow\_{SW \mid D}) \tag{8}$$

where *a* is the albedo.

The albedo depends strongly on the wave length, but generally a mean value is used for the whole visible spectrum. Besides its dependance on the surface type, it is important to consider its dependence on soil water content, vegetation state and surface roughness. The albedo depends moreover on the sun rays inclination, in particular for smooth surfaces: for small angles it increases. There is very rich literature about albedo description, it being a key parameter in the radiative exchange models, see for example Kondratyev (1969). Albedo is often divided in visible, near infrared, direct and diffuse albedo, as in Bonan (1996).

## **2.1.5 Long-wave radiation emitted by the surface** *R* ↑*LW*

This term indicates the long-wave radiation emitted by the earth surface, considered as a grey body with emissivity *ε<sup>s</sup>* (values from 0.95 to 0.98). The surface temperature *Ts* [*K*] is unknown and must be calculated by a LSM. *<sup>σ</sup>* <sup>=</sup> 5.6704 · <sup>10</sup>−<sup>8</sup> *<sup>W</sup>*/(*m*2*K*4) is Stefan-Boltzman constant.

$$R \uparrow\_{LW} = \varepsilon\_s \sigma T\_s^4 \tag{9}$$

## **2.1.6 Reflected long-wave radiation** *R* ↑*LW R*

This term is small and can be subtracted by the incoming long-wave radiation, assuming surface emissivity *ε<sup>s</sup>* equal to surface absorptivity:

$$R \downarrow\_{LW} = \varepsilon\_s \cdot \varepsilon\_a \sigma T\_a^4 \tag{10}$$

## **2.1.7 Radiation emitted and reflected by surrounding surfaces** *R* ↓*SW O* +*R* ↓*LW O*

It indicates the radiation reflected (*R* ↑*SW* +*R* ↑*LW R*) and emitted (*R* ↑*LW*) by the surfaces adjacent to the point considered. This term is important at small scale, in the presence of artificial obstructions or in the case of a very uneven orography. To calculate it with precision it is necessary to consider reciprocal orientation, illumination, emissivity and the albedo of every element, through a recurring procedure (Helbig et al., 2009). A simple solution is proposed for example in Bertoldi et al. (2005).

If the intervisible surfaces are hypothesized to be in radiative equilibrium, i.e. they absorb as much as they emit, these terms can be quantified in a simplified way:

$$\begin{array}{l} \mathbb{R}\downarrow\_{SW\ O} = (1-V)\mathbb{R}\uparrow\_{SW} \\ \mathbb{R}\downarrow\_{LW\ O} = (1-V)(\mathbb{R}\uparrow\_{LW} + \mathbb{R}\uparrow\_{LW\ R}) \end{array} \tag{11}$$

#### **2.1.8 Net radiation**

4 Will-be-set-by-IN-TECH

and there exist many models providing the soil incident radiation spectrum in a detailed way,

This term is a function of the atmospheric radiation due to Rayleigh dispersion and to the aereosols dispersion, as well as to the presence of cloud cover. The *R* ↓*SW D* actually is not isotropic and it depends on the sun position above the horizon. For its parametrization, see,

Often this term in not provided by standard meteorological measurements, and many LSMs provide expressions to calculate it. This term indicates the long-wave radiation emitted by atmosphere towards the earth. It can be calculated starting from the knowledge of the distribution of temperature, humidity and carbon dioxide of the air column above. If this information is not available, various formulas, based only on ground measurements, can be

*<sup>R</sup>* <sup>↓</sup>*LW D*<sup>=</sup> *�aσT*<sup>4</sup>

Usually for *�<sup>a</sup>* empirical formulas have been used, but it is also possible to provide a derivation based on physical topics like in Prata (1996). The cloud cover effect on this term is significant and not easy to consider in a simple way. Cloud cover data can be provided during the day

The albedo depends strongly on the wave length, but generally a mean value is used for the whole visible spectrum. Besides its dependance on the surface type, it is important to consider its dependence on soil water content, vegetation state and surface roughness. The albedo depends moreover on the sun rays inclination, in particular for smooth surfaces: for small angles it increases. There is very rich literature about albedo description, it being a key parameter in the radiative exchange models, see for example Kondratyev (1969). Albedo is

This term indicates the long-wave radiation emitted by the earth surface, considered as a grey body with emissivity *ε<sup>s</sup>* (values from 0.95 to 0.98). The surface temperature *Ts* [*K*] is unknown

often divided in visible, near infrared, direct and diffuse albedo, as in Bonan (1996).

**2.1.5 Long-wave radiation emitted by the surface** *R* ↑*LW*

by ground or satellite observations but, especially on night, is difficult to collect.

*<sup>a</sup>* (7)

*R* ↑*SW*= *a*(*R* ↓*SW P* +*R* ↓*SW D*) (8)

considering the various attenuation effects separately (Kondratyev, 1969).

**2.1.2 Diffuse downward short-wave radiation** *R* ↓*SW D*

**2.1.3 Diffuse downward long-wave radiation** *R* ↓*LW D*

*�<sup>a</sup>* atmosphere emissivity *f*(*ea*, *Ta*, cloud cover);

found in literature with expressions as follows:

**2.1.4 Reflected short-wave radiation** *R* ↑*SW* This term indicates the short-wave energy reflection.

for example, Paltrige & Platt (1976).

*Ta* air temperature [K];

*ea* vapor pressure [mb];

where *a* is the albedo.

with:

Inserting expressions (7) and ( 9) in the (3), the net radiation is:

$$R\_{\rm ll} = [sw \cdot R \downarrow\_{SW \, P} + V \cdot R \downarrow\_{SW \, D}](1 - V \cdot a) + V \cdot \varepsilon\_s \cdot \varepsilon\_{\rm d} \cdot \sigma \cdot T\_{\rm d}^4 - V \cdot \varepsilon\_s \cdot \sigma \cdot T\_s^4 \tag{12}$$

with *ε<sup>a</sup>* = *f*(*ea*, *Ta*, cloud cover) as for example in Brutsaert (1975).

Equation (12) is not invariant with respect to the spatial scale of integration: indeed it contains non-linear terms in *Ta*, *Ts*, *ea*, consequently the same results are not obtained if the local values of these quantities are substituted by the mean values of a certain surface. Therefore, the shift from a treatment valid at local level to a distributed model valid over a certain spatial scale must be done with a certain caution.

#### **2.1.9 Radiation adsorption and backscattering by vegetation**

Expression (12) needs to be modified to take into account the radiation adsorption and backscattering by vegetation, as shown in Figure 2. This effect is very important to obtain a correct soil surface skin temperature (Deardorff, 1978). From Best (1998) it is possible to derive the following relationship:

$$\begin{aligned} R\_{\mathcal{H}} &= \left[ \boldsymbol{sw} \cdot \boldsymbol{R} \downarrow\_{SW \cdot P} + \boldsymbol{V} \cdot \boldsymbol{R} \downarrow\_{SW \cdot D} \right] \left( \mathbf{1} - \boldsymbol{V} \cdot \boldsymbol{a} \right) \* \left( f\_{\text{trans}} + \boldsymbol{a}\_{\upsilon} \right) \\ &+ \left( \mathbf{1} - \boldsymbol{\varepsilon}\_{\upsilon} \right) \cdot \boldsymbol{V} \cdot \boldsymbol{\varepsilon}\_{\texttt{s}} \cdot \boldsymbol{\varepsilon}\_{\texttt{a}} \cdot \boldsymbol{\sigma} \cdot \boldsymbol{T}\_{\texttt{a}}^{\texttt{4}} + \boldsymbol{\varepsilon}\_{\upsilon} \cdot \boldsymbol{\varepsilon}\_{\texttt{s}} \cdot \boldsymbol{\sigma} \cdot \boldsymbol{T}\_{\upsilon}^{\texttt{4}} \end{aligned} \tag{13}$$

where *Tv* is vegetation temperature, *ε<sup>v</sup>* vegetation emissivity (supposed equal to absorption), *av* vegetation albedo (downward albedo supposed equal to upward albedo) and *ftrasm*

vegetation transmissivity, depending on plant type, leaf area index and photosynthetic activity.

Models oriented versus ecological applications have a very detailed parametrization of this term (Dickinson et al., 1986). Bonan (1996) uses a two-layers canopy model. Law et al. (1999) explicit the relationship between leaf area distribution and radiative transfer. A first energy budget is made at the canopy cover layer, and the energy fluxes are solved to find the canopy temperature, then a second energy budget is made at the soil surface. Usually a fraction of the grid cell is supposed covered by canopy and another fraction by bare ground.

Fig. 2. Schematic diagram of short-wave radiation (left) and long-wave radiation (right) absorbed, transmitted and reflected by vegetation and ground , as in equation 13 (from Bonan (1996), modified).

#### **2.2 Soil heat flux**

The soil heat flux *G* at a certain depth *z* depends on the temperature gradient as follows:

$$G = -\lambda\_s \frac{\partial T\_s}{\partial z} \tag{14}$$

where *λ<sup>s</sup>* is the soil thermal conductivity (*λ<sup>s</sup>* = *ρscsκ<sup>s</sup>* with *ρ<sup>s</sup>* density, *cs* specific heat and *κ<sup>s</sup>* soil thermal diffusivity) depending strongly on the soil saturation degree. The heat transfer inside the soil can be described in first approximation with Fourier conduction law:

$$\frac{\partial T\_s}{\partial t} = \kappa\_s \frac{\partial^2 T\_s}{\partial z^2} \tag{15}$$

Equation (14) neglects the heat associated to the vapor transportation due to a vertical gradient of the soil humidity content as well as the horizontal heat conduction in the soil. The vapor transportation can be important in the case of dry climates (Saravanapavan & Salvucci, 2000). The soil heat flux can be calculated with different degrees of complexity. The most simple assumption (common in weather forecast models) is to calculate *G* as a fraction of net radiation (Stull (1988) suggests *G* = 0.1*Rn*). Another simple approach is to use the analytical solution for a sinusoidal temperature wave. A compromise between precision and computational work is the force restore method (Deardorff, 1978; Montaldo & Albertson, 2001), still used in many hydrological models (Mengelkamp et al., 1999). The main advantage is that only two soil layers have to be defined: a surface thin layer, and a layer getting down to a depth where the daily flux is almost zero. The method uses some results of the analytical solution for a sinusoidal forcing and therefore, in the case of days with irregular temperature trend, it provides less precise results.

The most general solution is the finite difference integration of the soil heat equation in a multilayered soil model (Daamen & Simmonds, 1997). However, this method is computationally demanding and it requires short time steps to assure numerical stability, given the non-linearity and stationarity of the surface energy budget, which is the upper boundary condition of the equation.

#### **2.2.1 Snowmelt and freezing soil**

In mountain environments snow-melt and freezing soil should be solved at the same time as soil heat flux. A simple snow melt model is presented in Zanotti et al. (2004), which has a lumped approach, using as state variable the internal energy of the snow-pack and of the first layer of soil. Other models consider a multi-layer parametrization of the snowpack (e.g. Bartelt & Lehning, 2002; Endrizzi et al., 2006). Snow interception by canopy is described for example in Bonan (1996). A state of the art freezing soil modeling approach can be found in Dall'Amico (2010) and Dall'Amico et al. (2011).

#### **2.3 Turbulent fluxes**

6 Will-be-set-by-IN-TECH

vegetation transmissivity, depending on plant type, leaf area index and photosynthetic

Models oriented versus ecological applications have a very detailed parametrization of this term (Dickinson et al., 1986). Bonan (1996) uses a two-layers canopy model. Law et al. (1999) explicit the relationship between leaf area distribution and radiative transfer. A first energy budget is made at the canopy cover layer, and the energy fluxes are solved to find the canopy temperature, then a second energy budget is made at the soil surface. Usually a fraction of the

**Canopy** 

**4**

**Tv**

**R**↓**LW = (1-**ε**v )R**↓**LW +** ε**v** σ **Tv**

**R**↑**LW =** ε**g** σ **Tg**

Fig. 2. Schematic diagram of short-wave radiation (left) and long-wave radiation (right) absorbed, transmitted and reflected by vegetation and ground , as in equation 13 (from

The soil heat flux *G* at a certain depth *z* depends on the temperature gradient as follows:

*G* = −*λ<sup>s</sup>*

inside the soil can be described in first approximation with Fourier conduction law:

*∂Ts <sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>κ</sup><sup>s</sup>*

where *λ<sup>s</sup>* is the soil thermal conductivity (*λ<sup>s</sup>* = *ρscsκ<sup>s</sup>* with *ρ<sup>s</sup>* density, *cs* specific heat and *κ<sup>s</sup>* soil thermal diffusivity) depending strongly on the soil saturation degree. The heat transfer

*∂Ts*

*∂*2*Ts*

**(1-**ε**v )R**↑**LW +** ε**v** σ **Tv**

**4**

**Ts**

**4+ (1-**ε**g )R**↓**LW**

**Ground** 

*<sup>∂</sup><sup>z</sup>* (14)

*<sup>∂</sup>z*<sup>2</sup> (15)

grid cell is supposed covered by canopy and another fraction by bare ground.

**Shortwave Longwave** 

**<sup>R</sup>**↓**SW atm av R**↓**SW <sup>R</sup>**↓**LW atm**

**Veg ads**

**av R**<sup>↓</sup> **<sup>f</sup> SW trasm R**↓**SW**

**ag (ftrasm + av) R**↓**SW**

Bonan (1996), modified).

**2.2 Soil heat flux**

activity.

A modeling of the ground heat and vapor fluxes cannot leave out of consideration the schematization of the atmospheric boundary layer (ABL), meant as the lower part of atmosphere where the earth surface properties influence directly the characteristics of the motion, which is turbulent. For a review see Brutsaert (1982); Garratt (1992); Stull (1988).

A flux of a passive tracer *x* in a turbulent field (as for example heat and vapor close to the ground), averaged on a suitable time interval, is composed of three terms: the first indicates the transportation due to the mean motion *v*, the second the turbulent transportation *x*� **v**� , the third the molecular diffusion *k*.

$$\overline{\mathbf{F}} = \overline{\mathbf{x}}\,\overline{\mathbf{v}} + \overline{\overline{\mathbf{x'}}\,\overline{\mathbf{v'}}} - k\nabla\mathbf{x} \tag{16}$$

The fluxes parametrization used in LSMs usually only considers as significant the turbulent term only. The molecular flux is not negligible only in the few centimeters close the surface, and the horizontal homogeneity hypothesis makes negligible the convective term.

#### **2.3.1 The conservation equations**

The first approximation done by all hydrological and LSMs in dealing with turbulent fluxes is considering the Atmospheric Boundary Layer (ABL) as subject to a stationary, uniform motion, parallel to a plane surface.

This assumption can become limitative if the grid size becomes comparable to the vertical heterogeneity scale (for example for a grid of 10 m and a canopy height of 10 m). In this situation horizontal turbulent fluxes become relevant. A possible approach is the Large Eddy Simulation (Albertson et al., 2001).

If previous assumptions are made, then the conservation equations assume the form:

• Specific humidity conservation, failing moisture sources and phase transitions:

$$k\_{\upsilon} \frac{\partial^2 \overline{q}}{\partial z^2} - \frac{\partial}{\partial z} (\overline{w' q'}) = 0 \tag{17}$$

where:

*kv* is the vapor molecular diffusion coefficient [*m*2/*s*]

*q* = *mv mv*+*md* is the specific humidity [vapor mass out of humid air mass].

• Energy conservation:

$$k\_h \frac{\partial^2 \overline{\theta}}{\partial z^2} - \frac{\partial}{\partial z} (\overline{w' \theta'}) - \frac{1}{\rho c\_p} \frac{\partial H\_R}{\partial z} = 0 \tag{18}$$

where:


$$-\frac{1}{\rho}\frac{\partial \overline{p}}{\partial \mathbf{x}} + 2\omega \sin \phi \,\overline{v} + \nu \frac{\partial^2 \overline{u}}{\partial z^2} - \frac{\partial}{\partial z} (\overline{w'u'}) = 0\tag{19}$$

$$-\frac{1}{\rho}\frac{\partial \overline{p}}{\partial y} - 2\omega \sin \phi \,\overline{u} + \nu \frac{\partial^2 \overline{v}}{\partial z^2} - \frac{\partial}{\partial z}(\overline{w'v'}) = 0\tag{20}$$

where:


The vertical motion equation can be reduced to the hydrostatic equation:

$$\frac{\partial p}{\partial z} = -\rho g.\tag{21}$$

In a turbulent motion the molecular transportation terms of the momentum, heat and vapor quantity, respectively *ν*, *kh* and *kv*, are several orders of magnitude smaller than Reynolds fluxes and can be neglected.

#### **2.3.2 Wind, heat and vapor profile at the surface**

8 Will-be-set-by-IN-TECH

The first approximation done by all hydrological and LSMs in dealing with turbulent fluxes is considering the Atmospheric Boundary Layer (ABL) as subject to a stationary, uniform

This assumption can become limitative if the grid size becomes comparable to the vertical heterogeneity scale (for example for a grid of 10 m and a canopy height of 10 m). In this situation horizontal turbulent fluxes become relevant. A possible approach is the Large Eddy

(*w*�*q*�) = 0 (17)

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>0</sup> (18)

(*w*�*u*�) = 0 (19)

(*w*�*v*�) = 0 (20)

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>−</sup>*ρg*. (21)

If previous assumptions are made, then the conservation equations assume the form: • Specific humidity conservation, failing moisture sources and phase transitions:

*mv*+*md* is the specific humidity [vapor mass out of humid air mass].

(*w*�*θ*�) <sup>−</sup> <sup>1</sup>

• The horizontal mean motion equations are obtained from the momentum conservation by

*∂*2*u <sup>∂</sup>z*<sup>2</sup> <sup>−</sup> *<sup>∂</sup> ∂z*

*∂*2*v <sup>∂</sup>z*<sup>2</sup> <sup>−</sup> *<sup>∂</sup> ∂z*

*ρcp*

*∂HR*

*kv ∂*2*q <sup>∂</sup>z*<sup>2</sup> <sup>−</sup> *<sup>∂</sup> ∂z*

*kv* is the vapor molecular diffusion coefficient [*m*2/*s*]

simplifying Reynolds equations (Stull, 1988; Brutsaert, 1982 cap.3):

+ 2*ω sin φ v* + *ν*

*<sup>∂</sup><sup>y</sup>* <sup>−</sup> <sup>2</sup>*<sup>ω</sup> sin <sup>φ</sup> <sup>u</sup>* <sup>+</sup> *<sup>ν</sup>*

The vertical motion equation can be reduced to the hydrostatic equation:

*∂p*

In a turbulent motion the molecular transportation terms of the momentum, heat and vapor quantity, respectively *ν*, *kh* and *kv*, are several orders of magnitude smaller than Reynolds

*kh ∂*2*θ <sup>∂</sup>z*<sup>2</sup> <sup>−</sup> *<sup>∂</sup> ∂z*

*kh* is the thermal diffusivity [*m*2/*s*] *HR* is the radiative flux [*W*/*m*2] *θ* is the potential temperature [K] *ρ* is the air density [*kg*/*m*3] *w* is the vertical velocity [*m*/*s*] .

> − 1 *ρ ∂p ∂x*

−1 *ρ ∂p*

*ω* is the earth angular rotation velocity [rad/s]

*ν* is the kinematic viscosity [*m*2/*s*]

*φ* is the latitude [rad] .

fluxes and can be neglected.

**2.3.1 The conservation equations**

motion, parallel to a plane surface.

Simulation (Albertson et al., 2001).

where:

where:

where:

*q* = *mv*

• Energy conservation:

Inside the ABL we can consider, with a good approximation, that the decrease in the fluxes intensity is linear with elevation. This means that in the first meters of the air column the fluxes and the friction velocity *u*∗ can be considered constant. Considering the momentum flux constant with elevation implies that also the wind direction does not change with elevation (in the layer closest to the soil, where the geostrofic forcing is negligible). In this way the alignment with the mean motion allows the use of only one component for the velocity vector, and the problem of mean quantities on uniform terrain becomes essentially one-dimensional, as these become functions of the only elevation *z*.

In the first centimeters of air the energy transportation is dominated by the molecular diffusion. Close to the soil there can be very strong temperature gradients, for example during a hot summer day. Soil can warm up much more quickly than air. The air temperature diminishes very rapidly through a very thin layer called *micro layer*, where the molecular processes are dominant. The strong ground gradients support the molecular conduction, while the gradients in the remaining part of the surface layer drive the turbulent diffusion. In the remaining part of the surface layer the potential temperature diminishes slowly with elevation.

The effective turbulent flux in the interface sublayer is the sum of molecular and turbulent fluxes. At the surface, where there is no perceptible turbulent flux, the effective flux is equal to the molecular one, and above the first cm the molecular contribution is neglegible. According to Stull (1988), the turbulent flux measured at a standard height of 2 m provides a good approximation of the effective ground turbulent flux.

Fig. 3. (a) The effective turbulent flux in diurnal convective conditions can be different from zero on the surface. (b) The effective flux is the sum of the turbulent flux and the molecular flux (from Stull, 1988).

Applying the concept of effective turbulent flux, the molecular diffusion term can be neglected, while the hypothesis of uniform and stationary limit layer leads to neglect the convective terms due to the mean vertical motion and the horizontal flux. The vertical flux at the surface can then be reduced to the turbulent term only:

$$
\overline{\mathbf{F\_z}} = \overline{x'w'} \tag{22}
$$

In the case of the water vapor, equation (17) shows that, if there is no condensation, the flux is:

$$ET = \lambda \rho \overline{w'q'}\tag{23}$$

where *ET* is the evaporation quantity at the surface, *ρ* the air density and *λ* is the latent heat of vaporization.

Similarly, as to sensible heat, equation (18) shows that the heat flux at the surface *H* is:

$$H = \rho c\_p \overline{w'} \overline{\theta'} \tag{24}$$

where *cp* is the air specific heat at constant pressure.

The entity of the fluctuating terms *w*�*u*� , *w*�*θ*� and *w*�*q*� remains unknown if further hypotheses (called closing hypotheses) about the nature of the turbulent motion are not introduced. The closing model adopted by the LSMs is Bousinnesq model: it assumes that the fluctuating terms can be expressed as a function of the vertical gradients of the quantities considered (diffusive closure).

$$
\tau\_{\mathfrak{X}} = -\rho u' w' = \rho \mathcal{K}\_M \partial \overline{u} / \partial z \tag{25}
$$

$$H = -\rho c\_p \overline{w'} \theta' = -\rho c\_p K\_H \partial \overline{\theta} / \partial z \tag{26}$$

$$ET = -\lambda \rho w' q' = -\rho \mathbf{K}\_W \partial \overline{q} / \partial z \tag{27}$$

where *KM* is the turbulent viscosity, *KH* and *KW* [*m*2/*s*] are turbulent diffusivity. Moreover a logarithmic velocity profile in atmospheric neutrality conditions is assumed:

$$\frac{k\,\,\mu}{\mu\_{\ast 0}} = \ln(\frac{z}{z\_o}) \tag{28}$$

where *k* is the Von Karman constant, *z*<sup>0</sup> is the aerodynamic roughness, evaluated in first approximation as a function of the height of the obstacles as *z*0/*hc* � 0.1 (for more precise estimates see Stull (1988) p.379; Brutsaert (1982) ch.5; Garratt (1992) p.87). In the case of compact obstacles (e.g. thick forests), the profile can be thought of as starting at a height *d*0, and the height *z* can be substituted with a fictitious height *z* − *d*0.


Table 1. Values of aerodynamic roughness length *z*<sup>0</sup> for various natural surfaces (from Brutsaert, 1982).

Also the other quantities *θ* and *q* have an analogous distribution. Using as scale quantities *θ*∗<sup>0</sup> = −*w*�*θ*� 0/*u*∗<sup>0</sup> e *q*∗<sup>0</sup> = −*w*�*q*� <sup>0</sup>/*u*∗<sup>0</sup> and substituting them in the (25), the following integration is obtained:

10 Will-be-set-by-IN-TECH

In the case of the water vapor, equation (17) shows that, if there is no condensation, the flux

where *ET* is the evaporation quantity at the surface, *ρ* the air density and *λ* is the latent heat

(called closing hypotheses) about the nature of the turbulent motion are not introduced. The closing model adopted by the LSMs is Bousinnesq model: it assumes that the fluctuating terms can be expressed as a function of the vertical gradients of the quantities considered (diffusive

where *KM* is the turbulent viscosity, *KH* and *KW* [*m*2/*s*] are turbulent diffusivity. Moreover a

<sup>=</sup> ln( *<sup>z</sup> zo*

where *k* is the Von Karman constant, *z*<sup>0</sup> is the aerodynamic roughness, evaluated in first approximation as a function of the height of the obstacles as *z*0/*hc* � 0.1 (for more precise estimates see Stull (1988) p.379; Brutsaert (1982) ch.5; Garratt (1992) p.87). In the case of compact obstacles (e.g. thick forests), the profile can be thought of as starting at a height

> Surface type *z*<sup>0</sup> [cm] Large water surfaces 0.01-0.06 Grass, height 1 cm 0.1 Grass, height 10 cm 2.3 Grass, height 50 cm 5 Vegetation, height 1-2 m 20 Trees, height 10-15 m 40-70 Big towns 165

logarithmic velocity profile in atmospheric neutrality conditions is assumed:

*d*0, and the height *z* can be substituted with a fictitious height *z* − *d*0.

0/*u*∗<sup>0</sup> e *q*∗<sup>0</sup> = −*w*�*q*�

*k u u*∗0

Table 1. Values of aerodynamic roughness length *z*<sup>0</sup> for various natural surfaces (from

Also the other quantities *θ* and *q* have an analogous distribution. Using as scale quantities

Similarly, as to sensible heat, equation (18) shows that the heat flux at the surface *H* is:

where *cp* is the air specific heat at constant pressure.

The entity of the fluctuating terms *w*�*u*�

*ET* = *λρw*�*q*� (23)

*H* = *ρcpw*�*θ*� (24)

, *w*�*θ*� and *w*�*q*� remains unknown if further hypotheses

*τ<sup>x</sup>* = −*ρu*�*w*� = *ρKM∂u*/*∂z* (25) *H* = −*ρcpw*�*θ*� = −*ρcpKH∂θ*/*∂z* (26) *ET* = −*λρw*�*q*� = −*ρKW∂q*/*∂z* (27)

<sup>0</sup>/*u*∗<sup>0</sup> and substituting them in the (25), the following

) (28)

is:

of vaporization.

closure).

Brutsaert, 1982).

*θ*∗<sup>0</sup> = −*w*�*θ*�

$$\frac{k(\theta - \theta\_0)}{\theta\_{\ast 0}} = \ln(\frac{z}{z\_T}) \tag{29}$$

$$\frac{k(q-q\_0)}{q\_{\ast 0}} = \ln(\frac{z}{z\_q}).\tag{30}$$

The boundary condition chosen is *θ* = *θ*<sup>0</sup> in *z* = *zT* and *q* = *q*<sup>0</sup> in *z* = *zq*. The temperature *θ*<sup>0</sup> then is not the ground temperature, but that at the elevation *zT*. The roughness height *zT* is the height where temperature assumes the value necessary to extrapolate a logarithmic profile. Analogously, *zq* is the elevation where the vapor concentration assumes the value necessary to extrapolate a logarithmic profile.

Indeed, close to the soil (interface sublayer) the logarithmic profile is not valid and then, to estimate *zT* and *zq*, it would be necessary to study in a detailed way the dynamics of the heat and mass transfer from the soil to the first meters of air.

If we consider a real surface instead of a single point, the detail requested to reconstruct accurately the air motion in the upper soil meters is impossible to obtain. Then there is a practical problem of difficult solution: on the one hand, the energy transfer mechanisms from the soil to the atmosphere operate on spatial scales of few meters and even of few cm, on the other hand models generally work with a spatial resolution ranging from tens of m (as in the case of our approach) to tens of km (in the case of mesoscale models). Models often apply to local scale the same parametrizations used for mesoscale. Therefore a careful validation test, even for established theories, is always important.

Observations and theory (Brutsaert, 1982, p.121) show that *zT* and *zq* generally have the same order of magnitude, while the ratio *zT <sup>z</sup>*<sup>0</sup> is roughly included between <sup>1</sup> <sup>5</sup> <sup>−</sup> <sup>1</sup> 10 .

#### **2.3.3 The atmospheric stability**

In conditions different from neutrality, when thermal stratification allows the development of buoyancies, Monin & Obukhov (1954) similarity theory is used in LSMs. The similarity theory wants to include the effects of thermal stratification in the description of turbulent transportation. The stability degree is expressed as a function of Monin-Obukhov length, defined as:

$$L\_{MO} = -\frac{u\_{0\*}^3 \theta\_0}{kg\overline{w'}\theta'}\tag{31}$$

where *θ*<sup>0</sup> is the potential temperature at the surface.

Expressions of the stability functions can be found in many texts of Physics of the Atmosphere, for example Katul & Parlange (1992); Parlange et al. (1995). The most known formulation is to be found in Businger et al. (1971). Yet stability is often expressed as a function of bulk Richardson number *RiB* between two reference heights, expresses as:

$$Ri\_B = \frac{\text{g } z \,\Delta\theta}{\overline{\theta}u^2} \tag{32}$$

where Δ*θ* is the potential temperature difference between two reference heights, and *θ* is the mean potential temperature.

If *RiB* > 0 atmosphere is steady, if *RiB* < 0 atmosphere is unsteady. Differently from *LMO*, *RiB* is also a function of the dimensionless variables *z*/*z*<sup>0</sup> e *z*/*zT*. The use of *RiB* has the advantage that it does not require an iterative scheme.

Expressions of the stability functions as a function of *RiB* are provided by Louis (1979) and more recently by Kot & Song (1998). Many LSMs use empirical functions to modify the wind profile inside the canopy cover.

From the soil up to an elevation *hd* = *f*(*z*0), limit of the interface sublayer, the logarithmic universal profile and Reynolds analogy are no more valid. For smooth surfaces the interface sublayer coincides with the viscous sublayer and the molecular transport becomes important. For rough surfaces the profile depends on the distribution of the elements present, in a way which is not easy to parametrize. Particular experimental relations can be used up to elevation *hd*, to connect them up with the logarithmic profile (Garratt, 1992, p. 90 and Brutsaert, 1982, p. 88). These are expressions of non-easy practical application and they are still little tested.

#### **2.3.4 Latent and sensible heat fluxes**

As consequence of the theory explained in the previous paragraph, the turbulent latent and sensible fluxes *H* and *LE* can be expressed as:

$$H = \rho c\_p w' \theta' = \rho c\_p \mathbb{C}\_H u(\theta\_0 - \theta) \tag{33}$$

$$ET = \lambda \rho \overline{w'q'} = \lambda \rho \mathbb{C}\_E u (q\_0 - q)\_\prime \tag{34}$$

where *θ*<sup>0</sup> − *θ* and *q*<sup>0</sup> − *q* are the difference between surface and measurement height of potential temperature and specific humidity respectively. *CH* and *CE* are usually assumed to be equal and depending on the bulk Richardson number (or on Monin-Obukhov lenght):

$$\mathbb{C}\_{H} = \mathbb{C}\_{Hn} F\_{H}(\mathrm{Ri}\_{B}) \, \tag{35}$$

where *CHn* is the heat bulk coefficient for neutral conditions:

$$\mathbb{C}\_{Hn} = \mathbb{C}\_{En} = \frac{k^2}{[\ln(z/z\_0)][\ln(z\_d/z\_T)]} \tag{36}$$

derived on Eq. 29 and depending on the wind speed *u*, the measurement height *z*, the temperature (or moisture) measurement height *za*, the momentum roughness length *z*<sup>0</sup> and the heat roughness length *zT*.

A common approach is the 'electrical resistance analogy' (Bonan, 1996), where the atmospheric resistance is expressed as:

$$r\_{aH} = r\_{aE} = (\mathbb{C}\_H \,\, u)^{-1} \tag{37}$$

#### **3. Evapotranspiration processes**

In order to convert latent heat flux in evapotranspiration the energy conservation must be solved at the same time as water mass budget. In fact, there must be a sufficient water quantity available for evaporation. Moreover, vegetation plays a key role.

#### **3.1 Unsaturated soil evaporation**

12 Will-be-set-by-IN-TECH

If *RiB* > 0 atmosphere is steady, if *RiB* < 0 atmosphere is unsteady. Differently from *LMO*, *RiB* is also a function of the dimensionless variables *z*/*z*<sup>0</sup> e *z*/*zT*. The use of *RiB* has the advantage

Expressions of the stability functions as a function of *RiB* are provided by Louis (1979) and more recently by Kot & Song (1998). Many LSMs use empirical functions to modify the wind

From the soil up to an elevation *hd* = *f*(*z*0), limit of the interface sublayer, the logarithmic universal profile and Reynolds analogy are no more valid. For smooth surfaces the interface sublayer coincides with the viscous sublayer and the molecular transport becomes important. For rough surfaces the profile depends on the distribution of the elements present, in a way which is not easy to parametrize. Particular experimental relations can be used up to elevation *hd*, to connect them up with the logarithmic profile (Garratt, 1992, p. 90 and Brutsaert, 1982, p. 88). These are expressions of non-easy practical application and they are still little tested.

As consequence of the theory explained in the previous paragraph, the turbulent latent and

where *θ*<sup>0</sup> − *θ* and *q*<sup>0</sup> − *q* are the difference between surface and measurement height of potential temperature and specific humidity respectively. *CH* and *CE* are usually assumed to be equal and depending on the bulk Richardson number (or on Monin-Obukhov lenght):

derived on Eq. 29 and depending on the wind speed *u*, the measurement height *z*, the temperature (or moisture) measurement height *za*, the momentum roughness length *z*<sup>0</sup> and

A common approach is the 'electrical resistance analogy' (Bonan, 1996), where the

In order to convert latent heat flux in evapotranspiration the energy conservation must be solved at the same time as water mass budget. In fact, there must be a sufficient water quantity

*CHn* <sup>=</sup> *CEn* <sup>=</sup> *<sup>k</sup>*<sup>2</sup>

*H* = *ρcpw*�*θ*� = *ρcpCHu*(*θ*<sup>0</sup> − *θ*) (33)

*ET* = *λρw*�*q*� = *λρCEu*(*q*<sup>0</sup> − *q*), (34)

*CH* = *CHnFH*(*RiB*), (35)

*raH* = *raE* = (*CH u*)−<sup>1</sup> (37)

[ln(*z*/*z*0)][ln(*za*/*zT*)] (36)

that it does not require an iterative scheme.

profile inside the canopy cover.

**2.3.4 Latent and sensible heat fluxes**

the heat roughness length *zT*.

atmospheric resistance is expressed as:

**3. Evapotranspiration processes**

sensible fluxes *H* and *LE* can be expressed as:

where *CHn* is the heat bulk coefficient for neutral conditions:

available for evaporation. Moreover, vegetation plays a key role.

If the availability of water supply permits to reach the surface saturation level, then evaporation is potential *ET* = *EP* and then we have air saturation at the surface *q*(*Ts*) = *q*∗(*Ts*) (the superscript <sup>∗</sup> stands here for saturation). If the soil is unsaturated, *q*(*Ts*) �= *q*∗(*Ts*) and different approaches are possible to quantify the water content at the surface, in dependance of the water budget scheme adopted.

1. A first possibility is to introduce then the concept of surface resistance *rg* to consider the moisture reduction with respect to the saturation value. As it follows from equation (34):

$$ET = \lambda \rho \mathbb{C}\_E \mu (q\_0 - q) = \lambda \rho \frac{1}{r\_a} (q\_0 - q) = \lambda \rho \frac{1}{r\_a + r\_\mathcal{g}} (q\_0^\* - q) \tag{38}$$

2. As an alternative, we can define a soil-surface relative moisture

$$r\_h = q\_0 / q\_0^\* \tag{39}$$

and then the expression for evaporation becomes:

$$ET = \lambda \rho \frac{1}{r\_a} (r\_h \, q\_0^\* - q) \tag{40}$$

An expression of *rh* as a function of the potential *ψ<sup>s</sup>* [m] (work required to extract water from the soil against the capillarity forces) and of the ratio of the soil water content *η* to the saturation water content *η<sup>s</sup>* is given in Philip & Vries (1957):

$$r\_h = \exp\left(-\left(\mathbf{g} / R\_\upsilon T\_s\right) \psi\_\mathbf{s} \left(\eta / \eta\_\mathbf{s}\right)^{-b}\right) \tag{41}$$

where *Rv* = 461.53 [*J*/(*kg K*)] is the gas constant for water vapor, *Ts* is the soil temperature, *b* an empirical constant. Tables of these parameters for different soil types can be found in Clapp & Hornberger (1978).

Another more simple expression frequently applied in models to link the value *rh* with the soil water content *η* is provided by Noilhan & Planton (1989):

$$r\_h = \begin{cases} 0.5(1 - \cos(\frac{\eta}{\eta\_k}\pi)) \text{ se } \eta < \eta\_k\\ 1 & \text{if } \eta \ge \eta\_k \end{cases} \tag{42}$$

where *η* is the moisture content of a soil layer with thickness *d*1, and *η<sup>k</sup>* is a critical value depending on the saturation water content: *η<sup>k</sup>* � 0.75*ηs*.

3. A third possibility, very used in large-scale models, is that of expressing the potential/real evaporation ration through a simple coefficient:

$$ET = \propto EP = \propto \lambda \rho \frac{1}{r\_a} (q\_0^\* - q) \tag{43}$$

The value of *x* can be connected to the soil water content *η* through the expression (Deardorff, 1978) (see Figure 4):

$$\mathbf{x} = \min(1, \frac{\eta}{\eta\_k}) \tag{44}$$

Fig. 4. Dependence of *x* and *rh* on the soil water content *η* (Eq. 44-42)

#### **3.2 Transpiration**

Usually transpiration takes into account the canopy resistance *rc* to add to the atmospheric resistance *ra*:

$$ET = \text{x } EP = \propto \lambda \rho \frac{1}{r\_a + r\_c} (q\_0^\* - q) \tag{45}$$

The canopy resistance depends on plant type, leaf area index, solar radiation, vapor pressure deficit, temperature and water content in the root layer. There is a wide literature regarding such dependence, see for example Feddes et al. (1978); Wigmosta et al. (1994).

Canopy interception and evaporation from wet leaves are important processes modeled that should be modelled, according to Deardorff (1978). It is possible to distinguish two fundamental approaches: single-layer canopy models and multi-layer canopy models.

#### **Single-layer canopy models** (or "big leaf" models)

The vegetation resistance is entirely determined by stomal resistance and only one temperature value, representative of both vegetation and soil, is considered. Moreover a vegetation interception function can be defined so as to define when the foliage is wet or when the evaporation is controlled by stomal resistance.

#### **Multi-layer canopy models**

These are more complex models in which a soil temperature *Tg*, different from the foliage temperature *Tf* , is considered. Therefore, two pairs of equations of latent and sensible heat flux transfer, from the soil level to the foliage level, and from the latter to the free atmosphere, must be considered (Best, 1998). Moreover the equation for the net radiation calculation must consider the energy absorption and the radiation reflection by the vegetation layer.

Deardorff (1978) is the first author who presents a two-layer model with a linear interpolation between zones covered with vegetation and bare soil, to be inserted into atmosphere general circulation models. Over the last years many detailed models have been developed, above all with the purpose of evaluating the *CO*<sup>2</sup> fluxes between vegetation and atmosphere. Particularly complex is the case of scattered vegetation, 14 Will-be-set-by-IN-TECH

Usually transpiration takes into account the canopy resistance *rc* to add to the atmospheric

The canopy resistance depends on plant type, leaf area index, solar radiation, vapor pressure deficit, temperature and water content in the root layer. There is a wide literature regarding

Canopy interception and evaporation from wet leaves are important processes modeled that should be modelled, according to Deardorff (1978). It is possible to distinguish two fundamental approaches: single-layer canopy models and multi-layer canopy models.

The vegetation resistance is entirely determined by stomal resistance and only one temperature value, representative of both vegetation and soil, is considered. Moreover a vegetation interception function can be defined so as to define when the foliage is wet or

These are more complex models in which a soil temperature *Tg*, different from the foliage temperature *Tf* , is considered. Therefore, two pairs of equations of latent and sensible heat flux transfer, from the soil level to the foliage level, and from the latter to the free atmosphere, must be considered (Best, 1998). Moreover the equation for the net radiation calculation must consider the energy absorption and the radiation reflection by

Deardorff (1978) is the first author who presents a two-layer model with a linear interpolation between zones covered with vegetation and bare soil, to be inserted into atmosphere general circulation models. Over the last years many detailed models have been developed, above all with the purpose of evaluating the *CO*<sup>2</sup> fluxes between vegetation and atmosphere. Particularly complex is the case of scattered vegetation,

*ra* + *rc*

(*q*∗

<sup>0</sup> − *q*) (45)

*ET* <sup>=</sup> *x EP* <sup>=</sup> *<sup>x</sup> λρ* <sup>1</sup>

such dependence, see for example Feddes et al. (1978); Wigmosta et al. (1994).

**Single-layer canopy models** (or "big leaf" models)

**Multi-layer canopy models**

the vegetation layer.

when the evaporation is controlled by stomal resistance.

Fig. 4. Dependence of *x* and *rh* on the soil water content *η* (Eq. 44-42)

**3.2 Transpiration**

resistance *ra*:

where evaporation is due to a combination of soil/vegetation effects, which cannot be schematized as a single layer (Scanlon & Albertson, 2003).

Fig. 5. **Above**: scheme representing a single-layer vegetation model. Linked both with atmosphere (with resistance *ra*) and with the deep soil (through evapotranspiration with resistance *rs*), vegetation and soil surface layer are assumed to have the same temperature *T*<sup>0</sup> *<sup>f</sup>* . **Below**: scheme representing a multilayer vegetation model. Linked both with atmosphere (with resistances *rb* and *ra*), and with the deep soil (through evapotranspiration with resistance *rs*), as well as with the soil under the vegetation (*rd*), vegetation and soil surface layer are assumed to have different temperature *Tf* and *Tg*. *Pg* is the rainfall reaching the soil surface (from Garratt, 1992).

Given the many uncertainties regarding the forcing data and the components involved (soil, atmosphere), and the numerous simplifying hypotheses, the detail requested in a vegetation cover scheme is not yet clear.

A single-layer description of vegetation cover (big-leaf) and a two-level description of soil represent probably the minimum level of detail requested. In general, if the horizontal scale is far larger than the vegetation scale, a single-layer model is sufficient (Garratt, 1992, p. 242), as in the case of the general circulation atmospheric models or of mesoscale hydrologic models for large basins. These models determine evaporation as if the vegetation cover were but a partially humid plane at the atmosphere basis. In an approach of this kind surface resistance, friction length, albedo and vegetation interception must be specified. The surface resistance must include the dependence on solar radiation or on soil moisture, as transpiration decreases when humidity becomes smaller than the withering point (Jarvis & Morrison, 1981). For the soil, different coefficients depending on moisture are requested, together with a functional relation of evaporation to the soil moisture.

#### **4. Water in soils**

Real evaporation is coupled to the infiltration process occurring in the soil, and its physically-based estimate cannot leave the estimation of soil water content consideration. The most simple schemes to account water in soils used in LSMs single-layer and two-layer methods. The most general approach, which allows water transport for unsaturated stratified soil, is based on the integration of Richards (1931) equation, under different degrees of approximations.

#### **4.1 Single layer or bucket method**

In this method the whole soil layer is considered as a bucket and real evaporation *E*<sup>0</sup> is a fraction *x* of potential evaporation *Ep*, with *x* proportional to the saturation of the whole soil.

$$E\_0 = \varkappa E\_p \tag{46}$$

with *x* expressed by Eq. (44). The main problem of this method is that evaporation does not respond to short precipitation, leading to surface saturation but not to a saturation of the whole soil layer (Manabe, 1969).

#### **4.2 Two-layer or force restore method**

This method is analogous to the one developed to calculate the soil heat flux, but it requires calibration parameters which are unlikely to be known. With this method it is possible to consider the water quantity used by plants for transpiration, considering a water extraction by roots in the deepest soil layer (Deardorff, 1978).

#### **4.3 Multilayer methods and Richards equation**

Richards (1931) equation and Darcy-Buckingham law govern the unsaturated water transport in isobar and isothermal conditions:

$$\vec{q} = -\mathbb{K}\nabla \ (z + \psi)\tag{47}$$

$$\frac{\partial \psi}{\partial \eta} \nabla \cdot (K \nabla \psi) - \frac{\partial K}{\partial z} = \frac{\partial \psi}{\partial t} \tag{48}$$

where *q* = (*qx*, *qy*, *qz*) is the specific discharge, *K* is the hydraulic conductivity tensor, *z* is the upward vertical coordinate and *ψ* is the suction potential or matrix potential.

The determination of the suction potential allows also a more correct schematization of the plant transpiration and it lets us describe properly flow phenomena from the water table to the surface, necessary to the maintenance of evaporation from the soils.

Richards equation is, rightfully, an energy balance equation, even if this is not evident in the modes from which it has been derived. Then the solutions of the equation (48) must be searched by assigning the water retention curve which relates *ψ* with the soil water content *η* and an explicit relation of the hydraulic conductivity as a function of *ψ* (or *η*). Both relationships depend on the type of terrain and are variable in every point. K augments with *η*, until it reaches the maximum value *Ks* which is reached at saturation.

Although the integration of the Richards equation is the only physically based approach, it requires remarkable computational effort because of the non linearity of the water retention curve. It is difficult to find a representative water retention curve because of the high degree of spatial variability in soil properties (Cordano & Rigon, 2008).

## **4.4 Spatial variability in soil moisture and evapotranspiration**

Topography controls the catchment-scale soil moisture distribution (Beven & Freer, 2001) and therefore water availability for ET. Two methods most frequently used to incorporate sub-grid variability in soil moisture and runoff production SVATs models are the variable infiltration capacity approach (Wood, 1991) and the topographic index approach (Beven & Kirkby, 1979). They represent computationally efficient ways to represent hydrologic processes within the context of regional and global modeling. A review and a comparison of the two methods can be found in Warrach et al. (2002).

More detailed approaches need to track surface or subsurface flow within a catchment explicitly. Such approaches, which require to couple the ET model with a distributed hydrological model, are particularly useful in mountain regions, as presented in the next section.

## **5. Evapotranspiration in Alpine Regions**

16 Will-be-set-by-IN-TECH

soil, different coefficients depending on moisture are requested, together with a functional

Real evaporation is coupled to the infiltration process occurring in the soil, and its physically-based estimate cannot leave the estimation of soil water content consideration. The most simple schemes to account water in soils used in LSMs single-layer and two-layer methods. The most general approach, which allows water transport for unsaturated stratified soil, is based on the integration of Richards (1931) equation, under different degrees of

In this method the whole soil layer is considered as a bucket and real evaporation *E*<sup>0</sup> is a fraction *x* of potential evaporation *Ep*, with *x* proportional to the saturation of the whole soil.

with *x* expressed by Eq. (44). The main problem of this method is that evaporation does not respond to short precipitation, leading to surface saturation but not to a saturation of the

This method is analogous to the one developed to calculate the soil heat flux, but it requires calibration parameters which are unlikely to be known. With this method it is possible to consider the water quantity used by plants for transpiration, considering a water extraction

Richards (1931) equation and Darcy-Buckingham law govern the unsaturated water transport

where *q* = (*qx*, *qy*, *qz*) is the specific discharge, *K* is the hydraulic conductivity tensor, *z* is the

The determination of the suction potential allows also a more correct schematization of the plant transpiration and it lets us describe properly flow phenomena from the water table to

Richards equation is, rightfully, an energy balance equation, even if this is not evident in the modes from which it has been derived. Then the solutions of the equation (48) must be searched by assigning the water retention curve which relates *ψ* with the soil water content *η* and an explicit relation of the hydraulic conductivity as a function of *ψ* (or *η*). Both relationships depend on the type of terrain and are variable in every point. K augments with

*∂η* ∇ · (*K*∇*ψ*) <sup>−</sup> *<sup>∂</sup><sup>K</sup>*

*E*<sup>0</sup> = *xEp* (46)

*q* = −*K*∇ (*z* + *ψ*) (47)

*<sup>∂</sup><sup>t</sup>* (48)

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> *∂ψ*

relation of evaporation to the soil moisture.

**4. Water in soils**

approximations.

**4.1 Single layer or bucket method**

whole soil layer (Manabe, 1969).

**4.2 Two-layer or force restore method**

in isobar and isothermal conditions:

by roots in the deepest soil layer (Deardorff, 1978).

*∂ψ*

upward vertical coordinate and *ψ* is the suction potential or matrix potential.

the surface, necessary to the maintenance of evaporation from the soils.

*η*, until it reaches the maximum value *Ks* which is reached at saturation.

**4.3 Multilayer methods and Richards equation**

In alpine areas, evapotranspiration (ET) spatial distribution is controlled by the complex interplay of topography, incoming radiation and atmospheric processes, as well as soil moisture distribution, different land covers and vegetation types.


Spatially distributed hydrological and land surface models (e.g., Ivanov et al., 2004; Kunstmann & Stadler, 2005; Wigmosta et al., 1994) are able to describe land surface interactions in complex terrain, both in the temporal and spatial domains. In the next section we show an example of the simulation of the ET spatial distribution in an Alpine catchment simulated with the hydrological model GEOtop (Endrizzi & Marsh, 2010; Rigon et al., 2006).

## **6. Evapotranspiration in the GEOtop model**

The GEOtop model describes the energy and mass exchanges between soil, vegetation and atmosphere. It takes account of land cover, soil moisture and the implications of topography on solar radiation. The model is open-source, and the code can be freely obtained from the web site: http://www.geotop.org/. There, we provide a brief description of the 0.875 version of the model (Bertoldi et al., 2005), used in this example. For details of the most recent numerical implementation, see Endrizzi & Marsh (2010).

The model has been proved to simulate realistic values for the spatial and temporal dynamics of soil moisture, evapotranspiration, snow cover (Zanotti et al., 2004) and runoff production, depending on soil properties, land cover, land use intensity and catchment morphology (Bertoldi et al., 2010; 2006).

The model is able to simulate the following processes: (i) coupled soil vertical water and energy budgets, through the resolution of the heat and Richard's equations, with temperature and water pressure as prognostic variables (ii) surface energy balance in complex topography, including shadows, shortwave and longwave radiation, turbulent fluxes of sensible and latent heat, as well as considering the effects of vegetation as a boundary condition of the heat equation (iii) ponding, infiltration, exfiltration, root water extraction as a boundary condition of Richard's equation (iv) subsurface lateral flow, solved explicitly and considered as a source/sink term of the vertical Richard's equation (v) surface runoff by kinematic wave, and (vi) multi-layer glacier and snow cover, with a solution of snow water and energy balance fully integrated with soil.

The incoming direct shortwave radiation is computed for each grid cell according to the local solar incidence angle, including shadowing (Iqbal, 1983). It is also split into a direct and diffuse component according to atmospheric and cloud transmissivity (Erbs et al., 1982). The diffuse incoming shortwave and longwave radiation is adjusted according to the theory described in Par. 2.1. The soil column is discretized in several layers of different thicknesses. The heat and Richards' equations are written respectively as:

$$\mathcal{L}\_l(P)\frac{\partial T}{\partial t} - \frac{\partial}{\partial z}\left[K\_l(P)\frac{\partial T}{\partial z}\right] = 0\tag{49}$$

$$\mathcal{L}\_{\hbar}(P)\frac{\partial P}{\partial t} - \frac{\partial}{\partial z}\left[\mathcal{K}\_{\hbar}(P)\left(\frac{\partial P}{\partial z} + 1\right)\right] - q\_{\rm s} = 0\tag{50}$$

Where *T* is soil temperature, *P* the water pressure, *Ct* the thermal capacity, *Kt* the thermal conductivity, *Ch* the specific volumetric storativity, *Kh* the hydraulic conductivity, and *qs* the source term associated with lateral flow. The variables *Ct*, *Kt*, *Ch*, and *Kh* depend on water content, and, in turn, on water pressure, and are therefore a source of non-linearity. At the bottom of the soil column a boundary condition of zero fluxes has been imposed.

The boundary conditions at the surface are consistent with the infiltration and surface energy balance, and are given in terms of surface fluxes of water (*Qh*) and heat (*Qt*) at the surface, namely:

$$Q\_h = \min\left[p\_{net}, K\_{h1}\frac{(h - P\_1)}{dz/2} + K\_{h1}\right] - E(T\_1, P\_1) \tag{51}$$

$$Q\_l = SW\_{\rm in} - SW\_{\rm out} + L\mathcal{W}\_{\rm in} - L\mathcal{W}\_{\rm out}(T\_1) - H(T\_1) - LE(T\_1) \tag{52}$$

Where *pnet* is the net precipitation, *Kh*<sup>1</sup> and *P*<sup>1</sup> are the hydraulic conductivity and water pressure of the first layer, h is the pressure of ponding water, *dz* the thickness of the first layer, *T*<sup>1</sup> the temperature of the first layer. *E* is evapotranspiration (as water flux), *SWin* and *SWout* are the incoming and outgoing shortwave radiation, *LWin* and *LWout* the incoming 18 Will-be-set-by-IN-TECH

the web site: http://www.geotop.org/. There, we provide a brief description of the 0.875 version of the model (Bertoldi et al., 2005), used in this example. For details of the most recent

The model has been proved to simulate realistic values for the spatial and temporal dynamics of soil moisture, evapotranspiration, snow cover (Zanotti et al., 2004) and runoff production, depending on soil properties, land cover, land use intensity and catchment morphology

The model is able to simulate the following processes: (i) coupled soil vertical water and energy budgets, through the resolution of the heat and Richard's equations, with temperature and water pressure as prognostic variables (ii) surface energy balance in complex topography, including shadows, shortwave and longwave radiation, turbulent fluxes of sensible and latent heat, as well as considering the effects of vegetation as a boundary condition of the heat equation (iii) ponding, infiltration, exfiltration, root water extraction as a boundary condition of Richard's equation (iv) subsurface lateral flow, solved explicitly and considered as a source/sink term of the vertical Richard's equation (v) surface runoff by kinematic wave, and (vi) multi-layer glacier and snow cover, with a solution of snow water and energy balance

The incoming direct shortwave radiation is computed for each grid cell according to the local solar incidence angle, including shadowing (Iqbal, 1983). It is also split into a direct and diffuse component according to atmospheric and cloud transmissivity (Erbs et al., 1982). The diffuse incoming shortwave and longwave radiation is adjusted according to the theory described in Par. 2.1. The soil column is discretized in several layers of different thicknesses.

> *∂T ∂z*

*Qt* = *SWin* − *SWout* + *LWin* − *LWout*(*T*1) − *H*(*T*1) − *LE*(*T*1) (52)

 *∂P ∂z* + 1 

Where *T* is soil temperature, *P* the water pressure, *Ct* the thermal capacity, *Kt* the thermal conductivity, *Ch* the specific volumetric storativity, *Kh* the hydraulic conductivity, and *qs* the source term associated with lateral flow. The variables *Ct*, *Kt*, *Ch*, and *Kh* depend on water content, and, in turn, on water pressure, and are therefore a source of non-linearity. At the

The boundary conditions at the surface are consistent with the infiltration and surface energy balance, and are given in terms of surface fluxes of water (*Qh*) and heat (*Qt*) at the surface,

> (*h* − *P*<sup>1</sup> *dz*/2 <sup>+</sup> *Kh*<sup>1</sup>

Where *pnet* is the net precipitation, *Kh*<sup>1</sup> and *P*<sup>1</sup> are the hydraulic conductivity and water pressure of the first layer, h is the pressure of ponding water, *dz* the thickness of the first layer, *T*<sup>1</sup> the temperature of the first layer. *E* is evapotranspiration (as water flux), *SWin* and *SWout* are the incoming and outgoing shortwave radiation, *LWin* and *LWout* the incoming

bottom of the soil column a boundary condition of zero fluxes has been imposed.

*pnet*, *Kh*<sup>1</sup>

= 0 (49)

− *qs* = 0 (50)

− *E*(*T*1, *P*1) (51)

numerical implementation, see Endrizzi & Marsh (2010).

The heat and Richards' equations are written respectively as:

*Ch*(*P*)

*Qh* = *min*

*Ct*(*P*) *∂T <sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>∂</sup> ∂z Kt*(*P*)

*∂P <sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>∂</sup> ∂z Kh*(*P*)

(Bertoldi et al., 2010; 2006).

fully integrated with soil.

namely:

and outgoing longwave radiation, *H* the sensible heat flux and *LE* the latent heat flux. *H* and *LE* are calculated taking into consideration the effects of atmospheric stability (Monin & Obukhov, 1954).

*E* is partitioned by evaporation or sublimation from the soil or snow surface *EG*, transpiration from the vegetation *ETC*, evaporation of the precipitation intercepted by the vegetation *EVC*. Every cell has a fraction covered by vegetation and a fraction covered by bare soil. In the 0.875 version of the model, a one-level model of vegetation is employed, as in Garratt (1992) and in Mengelkamp et al. (1999): only one temperature is assumed to be representative of both soil and vegetation. In the most recent version, a two-layer canopy model has been introduced. Bare soil evaporation *Eg* is related to the water content of the first layer through the soil resistance analogy (Bonan, 1996):

$$E\_G = (1 - \alpha p) \ E\_P \, \frac{r\_d}{r\_d + r\_s} \tag{53}$$

where *cop* is fraction of soil covered by the vegetation *EP* is the potential *ET* calculated with equation 34 and *ra* the aerodynamic resistance:

$$r\_a = 1/\left(\rho \,\mathbb{C}\_E \,\mathbb{A}\right) \tag{54}$$

The soil resistance *rs* is function of the water content of the first layer.

$$r\_s = r\_a \frac{1.0 - (\eta\_1 - \eta\_r) / (\eta\_s - \eta\_r)}{(\eta\_1 - \eta\_r) / (\eta\_s - \eta\_r)}\tag{55}$$

where *η*<sup>1</sup> is the water content of the first soil layer close the surface, *η<sup>r</sup>* is the residual water content (defined following Van Genuchten, 1980) and *η<sup>s</sup>* is the saturated water content, both in the first soil layer.

The evaporation from wet vegetation is calculated following Deardorff (1978):

$$E\_{V\mathbb{C}} = \operatorname{cap} E\_P \,\delta\_W \tag{56}$$

where *δ<sup>W</sup>* is the wet vegetation fraction.

The transpiration from dry vegetation is calculated as:

$$E\_{T\mathbb{C}} = \left\langle \text{op } E\_P(1 - \delta\_W) \sum\_{i}^{n} \frac{f\_{\text{root}}^{i} r\_a}{r\_a + r\_c^i} \right\rangle \tag{57}$$

The root fraction *f <sup>i</sup> root* of each soil layer *i* is calculated decreasing linearly from the surface to a maximum root depth, depending from the cover type. The canopy resistance *rc* depends on solar radiation, vapor pressure deficit, temperature as in Best (1998) and on water content in the root zone as in Wigmosta et al. (1994).

#### **6.1 The energy balance at small basin scale: application to the Serraia Lake.**

An application of the model to a small basin is shown here, in order to bring out the problems arising when passing from local one-dimensional scale to basin-scale. The Serraia Lake basin is a mountain basin of 9 *km*2, with an elevation ranging from 900 to 1900 *m*, located in Trentino, Italy. Within the basin there is a lake of about 0.5 *km*2. During the year 2000 a study to calculate the yearly water balance was performed (Bertola & al., 2002).

The model was forced with meteorological measurement of a station located in the lower part of the basin at about 1000 m, and the stream-flow was calibrated for the sub-catchment of Foss Grand, of about 4 *km*2. Then the model was applied to the whole basin. Further details on the calibration can be found in Salvaterra (2001). Meteorological data are assumed to be constant across the basin, except for temperature, which varies linearly with elevation (0.6 *oC* / 100 m) and solar radiation, which slightly increases with elevation and is affected by shadow and aspect.

With the *GEOTOP* model it is possible to simulate the water and energy balance, aggregated for the whole basin (see figure 6 and 7) and its distribution across the basin. Figure 7 shows the map of the seasonal latent heat flux (ET) in the basin. During winter and fall ET is low (less than 40 *W*/*m*2), with the lowest values in drier convex areas. During summer and spring ET increases (up to 120 *W*/*m*2), with highest values in the bottom of the main valley (where indeed there are a lake and a wetland) and lowest values in north-facing, high-elevation areas.

Fig. 6. Monthly energy balance for the Serraia basin (TN, Italy).

The main factors controlling the ET pattern in a mountain environment (see Figure 8) are also: elevation, which controls temperature, aspect, which influences radiation, soil thickness, which determines storage capacity, topographic convergence, which controls the moisture availability. In particular, aspect has a primary effect on net radiation and a secondary effect more on sensible rather than on latent heat flux, as in Figure 9, where south aspect locations have larger *Rn* and *H*, but similar behavior for the other energy budget components). Water content changes essentially the rate between latent and heat flux, as in Figure 10 where wet locations have larger *ET* and lower *H*.

Therefore, the surface fluxes distribution seems to agree with experience and current hydrology theory, but the high degree of variability poses some relevant issues because the hypothesis of homogenous turbulence at the basis of the fluxes calculation is no more valid 20 Will-be-set-by-IN-TECH

The model was forced with meteorological measurement of a station located in the lower part of the basin at about 1000 m, and the stream-flow was calibrated for the sub-catchment of Foss Grand, of about 4 *km*2. Then the model was applied to the whole basin. Further details on the calibration can be found in Salvaterra (2001). Meteorological data are assumed to be constant across the basin, except for temperature, which varies linearly with elevation (0.6 *oC* / 100 m) and solar radiation, which slightly increases with elevation and is affected by shadow and

With the *GEOTOP* model it is possible to simulate the water and energy balance, aggregated for the whole basin (see figure 6 and 7) and its distribution across the basin. Figure 7 shows the map of the seasonal latent heat flux (ET) in the basin. During winter and fall ET is low (less than 40 *W*/*m*2), with the lowest values in drier convex areas. During summer and spring ET increases (up to 120 *W*/*m*2), with highest values in the bottom of the main valley (where indeed there are a lake and a wetland) and lowest values in north-facing, high-elevation areas.

**W/m2 Monthly basin averaged Rn <sup>G</sup> <sup>H</sup> ET**

**Flussi medi mensili per l' intero bacino**

**energy fluxes**

Fig. 6. Monthly energy balance for the Serraia basin (TN, Italy).

**-50**

locations have larger *ET* and lower *H*.

**0**

**50**

**100**

**150**

**200**

Energy balance: October 1999 - August 2000

The main factors controlling the ET pattern in a mountain environment (see Figure 8) are also: elevation, which controls temperature, aspect, which influences radiation, soil thickness, which determines storage capacity, topographic convergence, which controls the moisture availability. In particular, aspect has a primary effect on net radiation and a secondary effect more on sensible rather than on latent heat flux, as in Figure 9, where south aspect locations have larger *Rn* and *H*, but similar behavior for the other energy budget components). Water content changes essentially the rate between latent and heat flux, as in Figure 10 where wet

Therefore, the surface fluxes distribution seems to agree with experience and current hydrology theory, but the high degree of variability poses some relevant issues because the hypothesis of homogenous turbulence at the basis of the fluxes calculation is no more valid

**ott-99 nov-99 dic-99 gen-00 feb-00 mar-00 apr-00 mag-00 giu-00 lug-00 ago-00**

aspect.

Fig. 7. Seasonal latent heat maps *ET* [*W*/*m*2] for the Serraia basin (TN, Italy).

Fig. 8. Example of evapotranspiration *ET* for the Serraia basin, Italy. Notice the elevation effect (areas more elevated have less evaporation); the aspect effect (more evaporation in southern slopes, left part of the image); the topographic convergence effect on water availability (at the bottom of the valley).

Fig. 9. Difference in energy balance between locations with the same properties but different aspect. Dotted lines are for a south aspect location, while continuos lines are for a north aspect location. It can be noticed how south aspect locations have larger *Rn* and *H*, but similar behavior for the other fluxes.

Fig. 10. Difference in energy balance between locations with the same properties but different soil saturation. Dotted lines are for a dry location, while continuos lines are for wet location. It can be noticed how wet locations have larger *ET* and lower *H*, but similar behavior for the other fluxes. The time lag in *Rn* is due to differences in aspect.

(Albertson & Parlange, 1999). Moreover, horizontal differences in surface fluxes can start local air circulations, which can affect temperature and wind surface values with a feedback effect. How much such processes may affect the energy and water balance of the whole basin is easy to quantify, but *GEOTOP* can be a powerful tool to explore these issues.

## **7. Conclusion**

22 Will-be-set-by-IN-TECH

**W/m<sup>2</sup> Rn 1 G 1 H 1 ET 1 dEdT 1 Rn G H ET**

**Flussi**

similar behavior for the other fluxes.

**Flussi**

other fluxes. The time lag in *Rn* is due to differences in aspect.

Differences in aspect :Continuous =N, Dot =S

Fig. 9. Difference in energy balance between locations with the same properties but different aspect. Dotted lines are for a south aspect location, while continuos lines are for a north aspect location. It can be noticed how south aspect locations have larger *Rn* and *H*, but

**W/m<sup>2</sup> Rn G H ET dE/dt**

**Rn G H ET dE/dt**

**16/5/00 0.00 17/5/00 0.00 18/5/00 0.00 19/5/00 0.00 20/5/00 0.00**

Differences in saturation:Continuous sat=0.95, Dot=0.35

Fig. 10. Difference in energy balance between locations with the same properties but different soil saturation. Dotted lines are for a dry location, while continuos lines are for wet location. It can be noticed how wet locations have larger *ET* and lower *H*, but similar behavior for the

**16/5/00 0.00 17/5/00 0.00 18/5/00 0.00 19/5/00 0.00 20/5/00 0.00**

This chapter illustrates the components of the energy budget needed to model evapotranspiration (ET) and provides an extended review of the fundamental equations and parametrizations available in the hydrological and land surface models literature. In alpine areas, ET spatial distribution is controlled by the complex interplay of topography, incoming radiation and atmospheric processes, as well as soil moisture distribution, different land covers and vegetation types. An application of the distributed hydrological model GEOtop to a small basin is shown here, in order to bring out the problems arising when passing from local one-dimensional scale to basin-scale ET models.

## **8. References**


24 Will-be-set-by-IN-TECH

Brooks, P. D. & Vivoni, E. R. (2008a). Mountain ecohydrology: quantifying the role

Brooks, P. & Vivoni, E. R. (2008b). Mountain ecohydrology: quantifying the role of vegetation in the water balance of montane catchments., *Ecohydrology* 1: 187–192. Brutsaert, W. (1975). On a derivable formula for long-wave radiation from clear skies, *Water*

Brutsaert, W. (1982). *Evaporation into the Atmosphere: Theory, Hystory and Applications*, Kluver

Businger, J. A., Wyngaard, J. C., Izumi, Y. & Bradley, E. F. (1971). Flux profile relationships in

Chow, F. K., Weigel, A. P., Street, R. L., Rotach, M. W. & Xue, M. (2006). High-resolution

verification, and sensitivity experiments, *J. Appl. Met. and Clim.* pp. 63–86. Clapp, R. B. & Hornberger, G. M. (1978). Empirical equations for some hydraulic properties,

Cordano, E. & Rigon, R. (2008). A perturbative view on the subsurface

Daamen, C. C. & Simmonds, L. P. (1997). Soil, water, energy and transpiration, a numerical

Dall'Amico, M. (2010). *Coupled water and heat transfer in permafrost modeling*, PhD thesis,

Dall'Amico, M., Endrizzi, S., Gruber, S. & Rigon, R. (2011). A robust and energy-conserving model of freezing variably-saturated soil, *The Cryosphere* 5(2): 469–484.

Deardorff, J. W. (1978). Efficient prediction of ground surface temperature and moisture with inclusion of a layer of vegetation, *J. Geophys. Res.* 83(C4): 1889–1903. Dickinson, R. E., Heanderson-Sellers, A., Kennedy, P. J. & Wilson, M. (1986). Biosphere

Dubayah, A., Dozier, J. & Davis, F. W. (1990). Topographic distribution of clear-sky radiation

Endrizzi, S., Bertoldi, G., Neteler, M. & Rigon, R. (2006). Snow cover patterns and evolution at

Erbs, D. G., Klein, S. A. & Duffie, J. A. (1982). Estimation of the diffuse radiation fraction for hourly, daily and monthly average global radiation., *Sol. Energy* 28(4): 293–304. Feddes, R., Kowalik, P. & Zaradny, H. (1978). Simulation of field water use and crop yield,

*of the 63rd Eastern Snow Conference*, Newark, Delaware USA, pp. 195–209. Endrizzi, S. & Marsh, P. (2010). Observations and modeling of turbulent fluxes during melt

over the Konza Prairie, Kansas, *Water Resour. Res.* 26(4): 679–690.

*Simulation Monographs*, PUDOC, Wageningen, p. 188pp.

Trento. Available from http://eprints-phd.biblio.unitn.it/335/.

large-eddy simulations of flow in a steep alpine valley. Part I: methodology,

water pressure response at hillslope scale, *Water Resour. Res.*

model of water and energy fluxes in soil profiles and sparse canopies, *Technichal*

Institute of Civil and Environmental Engineering, Universita' degli Studi di Trento,

Atmosphere Transfer Scheme (BATS) for the NCAR Community Climate Model,

basin scale: GEOtop model simulations and remote sensing observations, *Proceedings*

at the shrub-tundra transition zone 1: point scale variations, *Hydrology Research*

the atmospheric surface layer, *J. Atmospheric Sciences* 28: 181–189.

10.1002/eco.27): 187 – 192.

*Resour. Res.* 11(5): 742–744.

*Water Resour. Res.* 14: 601–605.

*report*, University of Reading.

41(6): 471–490.

44(W05407): doi:10.1029/2006WR005740.

URL: *http://www.the-cryosphere.net/5/469/2011/*

*Technical Note NCAR/TN-275+STR*, NCAR.

Academic Publisher.

of vegetation in the water balance of montane catchments, *Ecohydrol.* 1(DOI:


Prata, A. J. (1996). A new long-wave formula for estimating downward clear-sky radiation at the surface., *Quarterly Journal of the Royal Meteorological Society* 122(doi: 10.1002/qj.49712253306): 1127–1151.

Richards, L. A. (1931). Capillary conduction of liquids in porous mediums, *Physics* 1: 318–333.

Rigon, R., Bertoldi, G. & Over, T. M. (2006). GEOtop: a distributed hydrological model with coupled water and energy budgets, *Journal of Hydrometeorology* 7: 371–388.

