Reference crop evapotranspiration estimate using high-resolution meteorological network ’ s data

Abstract. Water management authorities need detailed information about each component of the hydrological balance. This document presents a method to estimate the evapotranspiration rate, initialized in order to obtain the reference crop evapotranspiration rate (ET0). By using an Optimal Interpolation (OI) scheme, the hourly observations of several meteorological variables, measured by a high-resolution local meteorological network, are interpolated over a regular grid. The analysed meteorological fields, containing detailed meteorological information, enter a model for turbulent heat fluxes estimation based on Monin-Obukhov surface layer similarity theory. The obtained ET0 fields are then post-processed and disseminated to the users.


Introduction
Lombardia's public environmental agency (ARPA Lombardia) has undertaken a number of projects aimed at improving knownledge of the water budget in order to support public authorities involved in water management.This document presents a procedure to estimate ET 0 , defined by Allen et al. (1998) as the evapotranspiration rate for the hypothetical reference crop (short grass with an ample supply of water).
ET 0 is a climatic parameter computed using meteorological data only, without considering crop characteristics and soil factors.
The procedure's outputs are ET 0 hourly fields (unit mm h −1 ).The estimated evapotranspiration rate can be aggregated in time and space.Finally, the ET 0 estimates are published in a hydrological bulletin, which is disseminated to the users.
The source of meteorological information is ARPA Lombardia's high-resolution meteorological network.The meteorological variables used are hourly averaged values of temperature, relative humidity, wind components, solar global radiation and hourly cumulated precipitation values.The observing system is composed of about three hundred automated weather stations.The stations spatial distribution presents strong inhomogeneities, moreover the sensor equipment vary from station to station.As an example, Fig. 1 Correspondence to: C. Lussana (c.lussana@arpalombardia.it) shows the spatial distribution of thermometers and pyranometers.
The spatial domain of interest consists of the part of Po plain inside Lombardia's administrative boundaries (see fields in Figs. 4 and 5).Elevations range from 10 m (in the south-eastern part) to less than 300 m (in the northern part).However, all the information available is used in the ET 0 estimation procedure in order to reduce border effects.
Observations from the meteorological network are interpolated over a regular grid with 1.5 km of resolution.For all variables except solar radiation, a statistical interpolation scheme is applied.At present, the interpolation of the hourly averaged solar global radiation is performed with a simple inverse distance weighting method, Fig. 4 shows an example.An effort is made to prevent observations affected by gross errors from entering the interpolation procedure.
It is worthwhile remarking that the grid spatial resolution accounts for orographic details, but the interpolated fields can correctly reproduce phenomena resolved by the spatial resolution of the observational network, i.e. tenths of kilometers.
Once obtained the meteorological fields, a model for turbulent heat fluxes estimation is applied at each grid point.The model setup is such that the obtained latent heat flux coincides with the energetic content associated to ET 0 .
In this document, Sect. 2 describes the statistical interpolation scheme, Sect. 3 presents the procedure for turbulent fluxes estimation and Sect. 4 describes how the procedure is adapted to estimate ET 0 .
Published by Copernicus Publications.

Statistical interpolation
The statistical interpolation scheme is an implementation of the Optimal Interpolation (OI; Gandin, 1963).OI produces the best (in the sense of minimum analysis error variance), linear, unbiased estimate of the atmospheric field.The statistical interpolation scheme is applied to hourly averaged values of temperature, relative humidity, wind and to hourly cumulated precipitation values.OI filters out the unresolved spatial scales.OI schemes often use a model-derived first guess, but the implementation reported here uses a background field built by observations detrending, thus the presented OI implementations are model-independent.Figure 2 shows temperature and relative humidity fields.
The description of the OI implementation for temperature and relative humidity can be found in Uboldi et al. (2008) while in Lussana et al. (2009) the application to other variables is discussed and several test cases are reported.At present, with respect to that work, here the OI implementation for precipitation uses raingauges measurements only, without the integration of radar-derived precipitation esti- mates.The background field is set to zero everywhere and the OI parameters maintain the analysis field as close as possible to the observations: an example is shown in Fig. 3.The precipitation analysis plays an important role in the ET 0 estimation procedure, but only as a trigger: where the analyzed field of precipitation has a value greater than 1 mm h −1 the ET 0 value is not computed.In the future, the integration of raingauges with radar-derived precipitation estimates will allow a better definition of the portion of the domain not interested by precipitation.

Turbulent fluxes estimation
This Section presents a procedure to determine the surface moisture flux, using single-level meteorological data only.The moisture flux is intended as representative of an area rather than of a single roughness element.Due to the threedimensional nature of turbulent vortices, the assumption of horizontal representativeness also imposes a constraint on the height above the ground where fluxes must be evaluated.
In fact, from the general Atmospheric Boundary Layer (ABL) theory, the Surface Layer (SL) (about 10% of the ABL vertical extension) can be divided in Roughness Sub-Layer (RSL) and Inertial Sub-Layer (ISL) (Garratt, 1994).The RSL -the lower portion of the SL -contains the canopy layer.On the one hand, turbulence in the RSL is influenced by individual roughness elements and the turbulencerelated variables exhibit strong gradients in the vertical direction.On the other hand, turbulence in the ISL is influenced by the integrated effect of many roughness elements, thus the turbulence-related variables are representative of a larger area and the values of the turbulent fluxes are assumed to be constant with height.Most of Lombardia's automatic stations are located within the ISL (in practice, only urban stations lie within the RSL).The meteorological fields entering the procedure presented in the current session refer then to the ISL.
Based on the outlined structure of the SL, the evapotranspiration rate is defined as the net rate of passage of water vapor across a horizontal reference plane inside the ISL.
Turbulence is the most efficient mechanism to exchange momentum, energy and mass inside the ABL.In the framework of turbulence theory it is more convenient to deal with the energetic content associated to the moisture flux: the turbulent latent heat flux H e (W m −2 ).The evaporation rate E (mm h −1 ) can be obtained by using: www.adv-sci-res.net/3/113/2009/Adv.Sci.Res., 3,[113][114][115][116][117][118]2009 where λ is the latent heat of evaporation and ρ w is the density of water.
The main assumptions taken for this work are summarized in the following.
In order to make turbulence more tractable, without loosing its main features, it is customary to assume that, locally, the ergodic condition applies (horizontal homogeneity and stationarity for the averaged values of the meteorological variables).Furthermore, it is assumed that the wind vanish at the aerodynamic roughness length z 0 -the bottom RSL boundary -as prescribed by the no-slip condition.Especially in case of tall canopies, also the zero-plane displacement height d must be introduced.Finally, the presence of ice and liquid water is not considered, then phase transitions are not taken into account.
The two turbulent heat fluxes, latent, H e , and sensible, H 0 , are related to the meteorological variables using the so-called scaling parameters: Where ρ is the air density, c p is the specific heat at constant pressure; u * , T * and q * are the scaling parameters for momentum, temperature and specific humidity, respectively.These last three parameters are (vertically) constant through all the ISL.
In analogy with Zdunkowski and Bott (2003), the Monin-Obukhov Similarity Theory (MOST) for the SL allows computing the scaling parameters through the iterative procedure: Here u, T and q are the wind velocity, temperature and specific humidity, referred to the height above the ground z (in the ISL); k is the von Karman constant; L * is the Monin-Obukhov length scale; n is the iteration index.The iteration stops when a prescribed accuracy for the L * estimate is reached.The functions Ψ are the characteristic functions for momentum (M) and heat (H) (Beljaars and Holtslag, 1991).The unknown temperature and specific humidity T 0 and q 0 are formally assigned to z 0 , but they must be intended as two parameters characterizing the whole RSL.
In order to evaluate Eqs. ( 4)-( 7), T 0 and q 0 are needed.To obtain their values, a number of approximations must be made.The most relevant is the hypothesis of proportionality between ISL's turbulent fluxes and the difference between RSL and ISL values for the correspondent meteorological variables: ) The resistances r H and r V have the physical units of s m −1 .Both Eqs. ( 8) and ( 9) implicitly assume steady conditions.Equation ( 9) takes into account the evaporation originating from saturated surfaces in the RSL; the specific humidity difference in Eq. ( 9) can be rewritten as a linearized function of temperature: Where ∆≡∂q s /∂T is evaluated at a reference temperature between T 0 and T .Moreover, δq≡q s (T ) − q.
Referring to Monteith (1981), the resistance r H in Eq. ( 8) is estimated by the aerodynamic resistance, r H r a , governing the diffusion of energy and masses between RSL and ISL.r V in Eq. ( 9) is estimated as r V r a +r s , where the surface resistance r s is introduced to account for the evaporation and transpiration processes in the RSL.In this work the "big-leaf" model is used: the canopy is treated as a plane located at the lower boundary of the RSL and the latent heat flux is determined by the evaporation of liquid water from the "big-leaf".With these assumptions, Eq. ( 9) is replaced by: where δq 0 ≡q s (T 0 ) − q 0 .A further relation between the turbulent heat fluxes can be obtained from the overall surface energy balance.By means of the prognostic equation for entalphy applied to the earth surface, the energy balance can be written as (Zdunkowski and Bott, 2003): Where R n is the net radiation and G is the heat flux stored in the soil.The left-hand side defines the available energy.Equation ( 12) is strictly valid for atmosphere-surface interface in case of bare soil.Nevertheless, Eq. ( 12) is intended as the energy balance equation at the canopy top because the horizontal flux of energy due to advection and the rate of energy storage per unit area in the canopy layer are neglected.Therefore, the turbulent heat fluxes in Eq. ( 12) are assumed to be fluxes in the ISL and are interpreted as surface fluxes.
The net radiation is estimated using the Net All-wave Radiation Parameterization model (NARP; Offerle et al., 2003) while for G the expression G=βR n is used (β=0.1 for daytime and β=0.5 nighttime).
By making use of Eqs. ( 8)-( 12), the turbulent heat fluxes are computed as (deRooy and Holtslag, 1999): Where γ=c p /λ is the psychrometric constant.Equations ( 13) and ( 14) are used to overcome the problem of eliminating T 0 and q 0 inside the iterative procedure of Eqs. ( 4)-( 7).The iterative procedure is then rewritten as: 18) For the iterative procedure, the three parameters z 0 , d and r s need to be specified.Moreover, the albedo of the surface is required by the NARP model.

ET 0 estimation
The general method presented in Sect. 3 allows estimating the turbulent heat fluxes from single-level meteorological data.In order to obtain ET 0 , the parameters in the procedure must be properly initialized.In Allen et al. (1998) the reference crop is defined as grass with height h = 0.12m and albedo equals to 0.23.Furthermore, the relation for the surface roughness parameters are indicated as z 0 =0.123 h and d=2/3 h.With respect to the r s value, in a review of the FAO method for hourly period Allen et al. ( 2006) suggest using r s =50 s m −1 during daytime and r s =200 s m −1 during nighttime.The ET 0 value is not computed for grid points where the hourly precipitation rate is greater than 1mm h −1 because in case of rain the assumptions made in Sect. 3 are not justified.Figure 5 shows a field of daily cumulated ET 0 , obtained by summing up the hourly ET 0 estimates.

Conclusions
The presented procedure for ET 0 estimation relies on the availability of interpolated meteorological fields and combines these fields with a model for turbulent fluxes estimation near the ground, properly initialized.
The amount of information that contribute to the production of a ET 0 hourly field is remarkable, in consideration of the domain extension.In fact, about one thousand of meteorological observations are used every hour.The quality of the input meteorological fields is of crucial importance.The OI implementations for the meteorological variables provide reliable and detailed meteorological fields.However, the estimate of solar global radiation over the domain must be improved.
The outputs of the turbulent fluxes estimation model must be compared with experimental data.Nevertheless, preliminary qualitative evaluations are very encouraging.
In principle, the turbulent fluxes estimation procedure could be set up to obtain the real evapotranspiration rate for the actual crop.This would imply a more complex treatment of both the available geographical information and the vegetation behaviour, as a consequence several complementary parameters would need to be tuned.It would be thus possible to compare hourly estimates with the evapotranspiration estimates obtained using, for example, the approach of Allen et al. (2006), then possibly obtain some insights about crop coefficient values.
Edited by: E. Koch Reviewed by: two anonymous referees

Figure 1 .
Figure 1.Orography and station locations in the Lombardia domain.Triangles: thermometers.Circles: pyranometers.The bold black line is the administrative boundary.The inset shows the geographical location of Lombardia, longitude 8.5 to 11.5 • E, latitude 44.6 to 45.7 • N