Leuning et al. 2008

1 downloads 0 Views 5MB Size Report
itudes of the sites vary from 54°N to 35°S, while climatic regions range from subarctic, .... et al., 2008; Tsai and Philpot, 1998]. The MODIS Lai provides ...... Bleue), Russ K.Monson (Niwot Ridge), Henry L. Gholz and Timothy A. Martin (Mize) ...
Click Here

WATER RESOURCES RESEARCH, VOL. 44, W10419, doi:10.1029/2007WR006562, 2008

for

Full Article

A simple surface conductance model to estimate regional evaporation using MODIS leaf area index and the Penman-Monteith equation Ray Leuning,1 Y. Q. Zhang,2 Amelie Rajaud,1 Helen Cleugh,1 and Kevin Tu3 Received 1 October 2007; revised 19 May 2008; accepted 30 July 2008; published 25 October 2008.

[1] We introduce a simple biophysical model for surface conductance, Gs, for use with

remotely sensed leaf area index (Lai) data and the Penman-Monteith (PM) equation to calculate daily average evaporation, E, at kilometer spatial resolution. The model for Gs has six parameters that represent canopy physiological processes and soil evaporation: gsx, maximum stomatal conductance; Q50 and D50, the values of solar radiation and atmospheric humidity deficit when the stomatal conductance is half its maximum; kQ and kA, extinction coefficients for visible radiation and available energy; and f, the ratio of soil evaporation to the equilibrium rate corresponding to the energy absorbed at the soil surface. Model parameters were estimated using 2–3 years of data from 15 flux station sites covering a wide range of climate and vegetation types globally. The PM estimates of E are best when all six parameters in the Gs model are optimized at each site, but there is no significant reduction in model performance when Q50, D50, kQ, and kA are held constant across sites and gsx and f are optimized (linear regression of modeled mean daily evaporation versus measurements: slope = 0.83, intercept = 0.22 mm/d, R2 = 0.80, and N = 10623). The average systematic root-mean-square error in daytime mean evaporation was 0.27 mm/d (range 0.09–0.50 mm/d) for the 15 sites. The average unsystematic component was 0.48 mm/d (range 0.28–0.71 mm/d). The new model for Gs with two parameters yields better estimates of E than an earlier, simple model Gs = cL Lai, where cL is an optimized parameter. Our study confirms that the PM equation provides reliable estimates of evaporation rates from land surfaces at daily time scales and kilometer space scales when remotely sensed leaf area indices are incorporated into a simple biophysical model for surface conductance. Developing remote sensing techniques to measure the temporal and spatial variation in f is expected to enhance the utility of the model proposed in this paper. Citation: Leuning, R., Y. Q. Zhang, A. Rajaud, H. Cleugh, and K. Tu (2008), A simple surface conductance model to estimate regional evaporation using MODIS leaf area index and the Penman-Monteith equation, Water Resour. Res., 44, W10419, doi:10.1029/2007WR006562.

1. Introduction [2] Accurate estimates of evaporation are required to reduce uncertainties in constructing weekly to monthly water balances at catchment and regional scales. Accurate estimates of water yield (runoff) is needed by water resource managers responsible for urban and rural water supplies and knowledge of soil water availability is required for applications such as agricultural production. Improved estimates of evaporation from soils and transpiration from vegetation (combined symbol E) can constrain both of these important quantities because E is the largest term in the terrestrial water balance after precipitation. Evaporation is also of interest to 1

CSIRO Marine and Atmospheric Research, Canberra, ACT, Australia. CSIRO Land and Water, Canberra, ACT, Australia. 3 Department of Integrative Biology and Center for Stable Isotope Biogeochemistry, University of California, Berkeley, California, USA. 2

Copyright 2008 by the American Geophysical Union. 0043-1397/08/2007WR006562$09.00

meteorological agencies because energy partitioning at the Earth’s surface affects weather and climate. A major problem is that these diverse applications typically require knowledge of distribution of E across catchments and regions at daily time scales, whereas measurements are mostly made at a point, such as at a flux tower or a stream flow gauging station from which the spatial distribution of E must be inferred. A key research question is therefore how to spatialize information from points to areas; combining in situ and remotely sensed observations offer part of the solution to this challenge [Cleugh et al., 2007; Zhang et al., 2008]. [3] Remotely sensed data from satellites provides temporally and spatially continuous information on biophysical properties of the land surfaces such as land surface temperature [Wan et al., 2002], leaf area index [Myneni et al., 2002], vegetation indices [Huete et al., 2002] and vegetation type [Los et al., 2000]. Given the challenges and needs described above, there has been considerable effort by the scientific community to use such remotely sensed data to estimate spatially and temporally distributed evaporation.

W10419

1 of 17

W10419

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

Several classes of evaporation models use remotely sensed radiative surface temperature to estimate E, such as SEBAL [Bastiaanssen et al., 1998a, 1998b], SEBS [Su, 2002, 2005], NTDI [McVicar and Jupp, 2002], the resistance surface energy balance (RSEB) [Kalma and Jupp, 1990], the triangle method [Nemani and Running, 1989; Nishida et al., 2003; Gillies and Carlson, 1995], and the dual-source model developed by Norman et al. [1995] and Kustas and Norman [1999]. [4] Cleugh et al. [2007] found that multiseasonal time series of E estimated using MODIS measurements of the radiative surface temperature and the resistance surface energy budget approach compared poorly to eddy flux measurements of evaporation for two Australian ecosystems. Instead, they obtained satisfactory results when E was estimated using the Penman-Monteith (PM) equation [Monteith, 1964; Thom, 1975] with a simple model for surface conductance to evaporation given by Gs = cL Lai + Gs,min where Lai is the leaf area index obtained from MODIS remote sensing and cL is a constant and Gs,min is the surface conductance controlling soil evaporation. They then used the PM algorithm to produce an evaporation climatology for Australia at 1 km resolution. Mu et al. [2007] also used the PM equation to estimate E, but found that the simple surface conductance model of Cleugh et al. [2007] was inadequate when tested against evaporation measurements from 19 AmeriFlux sites. Mu et al. [2007] revised the model for Gs by introducing scaling functions that range between 0 and 1 to account for the response of stomata to humidity deficit of the air, Da, and air temperature, Ta. They also introduced a separate term for evaporation from the soil surface, a term that was acknowledged by Cleugh et al. [2007], but which they incorporated via the coefficient cL. The revised Gs algorithm of Mu et al. [2007] resulted in good agreement between predictions of E by the PM equation and the AmeriFlux measurements. They then used the revised algorithm to predict the global distribution of mean annual E at a spatial resolution of 0.05°. [5] The results of Cleugh et al. [2007] and Mu et al. [2007] show that the PM equation is a biophysically sound and robust framework for estimating daily E at regional to global scales using remotely sensed data. This paper builds on the original proposition of Cleugh et al. [2007] and the important developments of Mu et al. [2007] by introducing a new, simple biophysical algorithm for Gs based on our understanding of leaf- and canopy-level plant physiology, radiation absorption by plant canopies and evaporation from the underlying soil surface. The formulation for Gs modifies that of Kelliher et al. [1995] (hereinafter referred to as K95), by including the influence of Da on stomatal conductance and by introducing a soil evaporation algorithm that is simpler than that proposed by Mu et al. [2007] and which has the potential to be estimated using remote sensing. Performance of the new algorithm for Gs is then tested on several years of evaporation measurements at each of 15 flux station sites covering a wide range in climate and vegetation globally, including deciduous and evergreen forests, a corn crop, a wetland, grassland and two woody savannas. In a companion paper, Zhang et al. [2008] use mean annual E estimated from runoff measured at the gauged catchments and leaf area indices obtained from MODIS remote sensing to estimate parameters of the Gs model developed in this

W10419

paper. They then used the calibrated model and the PM equation to provide spatially resolved estimates of E across the Murray Darling Basin in SE Australia. [6] The following section describes the surface conductance model, followed by a brief description of the data sources and methods used to control the quality of the data. We then examine the performance of the model as we reduce from six to two the number of parameters to be estimated for each vegetation type. Model performance is assessed against measurements at the individual flux stations and across all vegetation types.

2. Theory 2.1. Surface Conductance [7] Since its derivation by Monteith [1964], the PenmanMonteith (PM) equation has been used extensively in the literature to describe the biophysics of evaporation, E, from land surfaces. The equation states that the flux of latent heat associated with E is given by   eA þ rcp =g Da Ga lE ¼ ; e þ 1 þ Ga =Gs

ð1Þ

where l is the latent heat of evaporation, e = s/g, in which g is the psychrometric constant and s = de*/dT is the slope of the curve relating saturation water vapor pressure to temperature, A is the available energy absorbed by the surface (net absorbed radiation, Rn minus soil heat flux, G), r is the density of air, cp is the specific heat of air at constant pressure, Da = e* (Ta)  ea is the water vapor pressure deficit of the air (humidity deficit), in which e* (Ta) is the saturation water vapor pressure at air temperature and ea is the actual water vapor pressure, Ga is the aerodynamic conductance and Gs is the surface conductance accounting for evaporation from the soil and transpiration by the vegetation. In this paper we assume that the quantities A, Da and Ga are either known or are estimated using measured meteorological fields of incoming solar radiation, temperature, humidity and wind speed [Cleugh et al., 2007]. The main challenge in practical application of equation (1) is thus to determine Gs. [8] Total evaporation is the sum of transpiration from the plant canopy, Ec, and evaporation from the soil, Es: E ¼ Ec þ Es :

ð2Þ

Using separate versions of the PM equation for E and Ec, and the assumption that evaporation from the soil occurs at some fraction, f, of the equilibrium rate at the soil surface, eAs/(1 + e), we can write     eA þ rcp =g Da Ga eAc þ rcp =g Da Ga f eAs : ¼ þ e þ 1 þ Ga =Gs e þ 1 þ Ga =Gc eþ1

ð3Þ

Note that we distinguish between surface conductance, Gs, and the canopy conductance Gc, and that we partition the total available energy A into that absorbed by the canopy, Ac, and by the soil, As. For the purposes of this study, f is a parameter to be estimated from the data, although in reality f will vary with moisture content of the soil near to the surface. Discussion of potential remote sensing techniques

2 of 17

W10419

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

W10419

Figure 1. Response of Gs/gsx to (a) Da, (b) A, (c) D50, (d) Q50, (e) gsx, and (f) f. Except when varied, parameter values are gsx = 0.008 m/s, Q50 = 30 W/m2, Ga = 0.033 m/s, kQ = kA = 0.6, Da = D50 = 1.0 kPa, f = 0.5, and A = 500 W/m2. Photosynthetically active radiation, Qh, is assumed to correlate with A according to Qh = 0.8A in all plots. for estimating the temporal and spatial variation of f is deferred until later. [9] The fraction of the total available energy absorbed by the canopy and by the soil are given respectively by Ac/A = 1  t and As/A = t, where t = exp (kALai), kA is an extinction coefficient for available energy and Lai is leaf area index. Note that this ignores the soil heat flux, G, because G ! 0 when we later use daily average meteorological variables. [10] With these definitions, equation (3) can be written as e þ Ga =Gi eð1  t Þ þ Ga =Gi f et ¼ þ ; e þ 1 þ Ga =Gs e þ 1 þ Ga =Gc eþ1

ð4Þ

Several useful limits can be derived from this equation. When f = 1, evaporation at the soil surface occurs at the equilibrium rate and equation (6) reduces to 2 3 tGa Ga 1 þ þ 6 ðe þ 1ÞGc eGi 7 7; Gs ¼ Gc 6 4 5 Ga 1tþ eGi

which is the same as equation (10) derived by K95 in their survey of maximum surface conductances of various land cover types. [11] When the soil surface is completely dry, f = 0 and then 2

3 Ga 6 7 eGi 7: Gs ¼ Gc 6 4 ðe þ 1ÞGc Ga 5 1þt þ Ga eGi

where Gi is the ‘‘climatological’’ conductance defined by Monteith [1964] as A  : Gi ¼  rcp =g Da

1tþ

ð5Þ

An expression for Gs in terms of Gc, t, f, Gi and Ga is obtained by rearranging equation (4):   2 3 tGa ðe þ 1Þð1  f ÞGc Ga þ 1 þ f  6 ðe þ 1ÞGc Ga eGi 7 7: ð6Þ   Gs ¼ Gc 6 4 5 ðe þ 1Þð1  f ÞGc Ga þ 1t f  Ga eGi

ð7Þ

ð8Þ

Note that Gs 6¼ Gc because Gc depends on radiation absorbed by the plant canopy, whereas Gs is a function of radiation absorbed by both canopy and soil. [12] If all radiation is absorbed by the canopy, t = 0, and thus

3 of 17

3 Ga 1 þ 6 eGi 7 7 ¼ Gc : GS ¼ Gc 6 4 Ga 5 1þ eGi 2

ð9Þ

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

W10419

W10419

Figure 2. Contour plots of the evaporative fraction, fE = lE/A, as a function of absorbed photosynthetically active radiation, Qh, and leaf area index (Lai 100) for (a) Da = 0.5 kPa, (b) Da = 1.0 kPa, (c) Da = 1.5 kPa, and (d) Da = 2.0 kPa. Fixed parameter values are gsx = 0.008 m/s, Q50 = 30 W/m2, Ga = 0.033 m/s, kQ = kA = 0.6, D50 = 1.0 kPa, and f = 0.5. Finally, when t = 1, Gc = 0 because there is no canopy, and then from equation (4) we obtain Gs ¼

efGa : ðe þ 1Þ½eð1  f Þ þ Ga =Gi

ð10Þ

In this case Ga and Gi appropriate to bare soil should be used, rather than those relevant to plant canopies. 2.2. Canopy Conductance and Evaporative Fraction [13] Gs is dependent on the as yet unspecified Gc in each of equations (6) – (9). K95 developed an expression for canopy conductance in terms of maximum stomatal conductance of leaves at the top of the canopy, gsx, Lai and an hyperbolic response to absorbed shortwave radiation. They showed that " # gsx Qh þ Q50   Gc ¼ ln ; kQ Qh exp kQ L þ Q50

ð11Þ

where Qh is the flux density of visible radiation at the top of the canopy (approximately half of incoming solar radiation), kQ is the extinction coefficient for shortwave radiation and Q50 is the visible radiation flux when stomatal conductance is half its maximum value. A similar expression was derived by Saugier and Katerji [1991] and Dolman et al. [1991]. Here we modify this expression to include the response of stomatal conductance to humidity deficit as developed by Leuning [1995] to give " #  gsx Qh þ Q50 1   ; Gc ¼ ln kQ Qh exp kQ L þ Q50 1 þ Da =D50

ð12Þ

where D50 is the humidity deficit at which stomatal conductance is half its maximum value. This formula was used by Isaac et al. [2004] and Wang et al. [2004], and found to be an improvement over that of K95. Substitution of equation (12) into (6) shows that model for Gs contains the six parameters gsx, Q50, D50 kQ, kA and f.

4 of 17

5 of 17

Tonzi

Tumbarumba

Virginia Park

1078

43

45

wsa

ebf

wsa

enf

Evergreen broadleaf Woody savanna

Slash pine clearcut Deciduous broadleaf forest Evergreen needleleaf forest Evergreen broadleaf forest Evergreen needleleaf forest Woody savanna

Wetland

Grassland

Deciduous broadleaf forest Deciduous broadleaf forest Mixed forest

Needleleaf forest

Croplands

IGBP Class

Quercus douglasii plus C3 winter grasses Eucalyptus delagatensis Eucalyptus creba, Eucalyptus, drepanophylla plus C4 grasses in summer

Picea mariana

Abies lasiocarpa, Picea engelmannii, Pinus contorta Cleared forest

Acer saccharum

Pinus elliottii

Sphagnum mosses

Red spruce, Eastern hemlock C4 grasses

Fagus sylvatica, Quercus petraea

Picea sitkensis, Pseudotsuga menziesii, Betula pendula Fagus sylvatica, Acer, Fraxinus

Corn/soybean

Dominant Species

0.5 –

1531 70

40 5–8

117 1200 200

38°25053.7600N, 120°57057.5400W 35°39020.600S, 148°907.500E 19°5300000S, 146°3301400E

USA

Australia

Australia

10.1

1 – 21

629

53°590 13.8100N, 105°70 4.0400W

Brazil Canada

11.5

27

275 3050

5

50



40°01058.3600N, 105°32047.0500W

19.5

60

45°12014.6500N, 68°4402500W 31°44011.5000N, 109°56030.7700W 45°24033.8400N, 75°3101200W 29°45053.2800N, 82°14041.3400W 39°19023.3400N, 86°24047.3000W

13

300

48°40027.1800N, 7°3052.6200E

33

430

9

340

51°4045.3600N, 10°2707.200E

3

300

40°0021.9600N, 88°17030.7200W 56°36023.5900N, 3°47048.5500W

3°104.90500S, 54°58017.16600W

USA

USA

USA

Canada

USA

USA

France

Germany

UK

USA

Country

4.0

4.7

2.6

39.8

19.9

31.5

38.9

28.9

42.9

35.2

36.9

26.6

34.8

15.7

36.2

30.7

37.6

38.1

17.6

9.6

33.7

32.6

23.6

14.0

30.8

22.8

5.5

13.2

34.0

Average

16.9

Minimum

21.8

9.2

16.0

2.2

26.0

1.1

13.0

20.0

6.8

13.1

7.2

10.5

8.0

7.2

11.7

Maximum

Temperature (°C)

402.9

1011.2

532.7

376.3

1598.5

562.9

1121.4

898.3

725.0

356

645.1

1154.5

853.9

1169.0

644.6

Annual Precipitation (mm)

References

Leuning et al. [2005]

Leuning et al. [2005]

Baldocchi et al. [2004]

Griffis et al. [2003], Mahrt and Vickers [2002]

Goulden et al. [2004], Miller et al. [2004]

Turnipseed et al. [2002]

Ehman et al. [2002], Oliphant et al. [2004]

Gholz and Clark [2002]

Admiral et al. [2006]



Hollinger et al. [1999]

Granier et al. [2000]

Knohl et al. [2003]

Falge et al. [2002]

Meyers and Hollinger [2004]

Full details on each site can be found at the ORNL DAAC by substituting the site identification (ID) number for xxx at http://www.fluxnet.ornl.gov/fluxnet/sitepage.cfm?SITEID=xxx.

SSA-Old Black Spruce

259

ebf

enf

dbf

bare

wet

grass

mf

dbf

dbf

enf

crp

Code

Canopy Height (m)

Elevation (m)

Latitude and Longitude

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

a

Santarem Km83

85

Mer Bleue

235

Niwot Ridge

Kendall

814

997

Howland

890

Morgan Munroe

Hesse

476

968

Hainich

392

Mize

Griffin

777

1071

Bondville

Site Name

830

Site ID

a

Table 1. Details of Fluxnet Sites Used to Evaluate Performance of the Surface Conductance Model in Estimating Land Surface Evaporation

W10419 W10419

W10419

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

W10419

Figure 3. (a) Mean annual temperature and precipitation at each of the 15 flux station sites used in the data analysis. (b) Mean Lai at each flux station versus mean annual precipitation. [14] Predictions of the ratio Gs/gsx, using equations (6) and (12), are shown as a function of Lai in Figure 1 for typical values of the meteorological parameters Da and A, the physiological parameters D50, Q50 and gsx and the soil evaporation parameter f. We note first that values of Gs/gsx are lower than in the work by K95 for all the sensitivity analyses in Figure 1 because the extra term 1/(1 + Da/D50) in equation (11) reduces Gc compared to corresponding values of K95. As a result, increasing Da from 0.5 to 2.0 kPa causes Gs/gsx to decrease by >50% at all Lai (Figure 1a), in contrast to the model of K95 that has no such dependence. Compared to the work by K95, the revised model also shows a greater dependence of Gs/gsx on A at high Lai (Figure 1b) because here we maintain a realistic correlation between Qh and A using Qh = 0.8A (R. Leuning, unpublished data, 2008), whereas K95 held Qh fixed when A was varied. Increasing values of the physiological parameter D50 causes Gs/gsx to increase, while increasing Q50 causes it to decrease at all Lai (Figures 1c and 1d). In contrast, the influence of varying gsx and f on Gs/gsx is only significant when Lai < 3 (Figures 1e and 1f). Soil evaporation is assumed to occur at the equilibrium rate in the K95 model (equivalent to f = 1) and this causes a minimum in Gs/gsx around Lai 1, whereas this minimum is absent or less pronounced in the new model, where f 1. Low rates of evaporation from sparse canopies correspond to low values of f and hence Gs. [15] Contour plots of the evaporative fraction, fE = lE/A, as predicted by the evaporation model, are shown in Figure 2 as a function of Qh and Lai 100 for values of Da = 0.5, 1.0, 1.5, and 2.0 kPa. Increasing Da for any combination of Qh and Lai results in an increase in evaporative fraction; for example, fE is predicted to increase from 0.75 to 1.45 for Qh = 300 W m2 and Lai = 3 as Da increases from 0.5 to 2.0 kPa. The contour spacing decreases with increasing Da, indicating increasing sensitivity of fE to variations in humidity deficit and leaf area index at higher humidity deficits. Sensitivity of fE to Lai is greatest at low values of Qh and vice versa.

[16] We next describe the data sources for evaporation fluxes, meteorology and Lai used to estimate model parameters and to evaluate model performance.

3. Methods 3.1. Sites, Fluxes, and Meteorological Data [17] Table 1 provides details of the 15 Fluxnet sites used to evaluate performance of the surface conductance model in estimating land surface evaporation. References describing each site are also listed. The sites are located across several continents and include savannas, grasslands, a corn crop, mixed deciduous forests, evergreen needleleaf and broadleaf forests, a wetland and cleared forest land. Latitudes of the sites vary from 54°N to 35°S, while climatic regions range from subarctic, cool temperate, warm temperate to tropical. Maritime and continental climates are also represented. Minimum temperatures at the continental SSA Old Black Spruce site reach 39.8°C while maximum temperatures exceeded 35°C at many sites. Annual precipitation across sites ranged from 376 to 1599 mm. [18] A plot of mean annual precipitation versus mean annual temperature (Figure 3a) shows that the sites chosen for analysis cover much of the observed range in global climate. The sites ranged from the boreal SSA-Old Black Spruce site with mean annual temperature of 1.2°C and precipitation of 376 mm to the tropical Santarem Km83 site with 1600 mm rainfall and mean annual temperatures of 26.1°C. High mean temperatures of 22.3°C are also observed at Virginia Park, which is a seasonally wet/dry savanna but with an annual rainfall of only 403 mm that occurs during the summer wet season. [19] For the analysis, time series of latent heat flux, wind speed, minimum and maximum air temperature, water vapor pressure, solar radiation, air pressure, net radiation and soil heat flux for the selected sites were obtained from the Oak Ridge National Laboratory, Distributed Active Archive Center (ORNL-DAAC, http://www.fluxnet.ornl. gov/MODIS/modis.html). [20] Latent heat flux data were eliminated from the analysis when measured values were consistently at zero

6 of 17

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

W10419

W10419

Figure 4. Scatterplots of measured daily average H + lE versus daily average Rn  G (=A) for each of the 15 flux sites used in this analysis. Slopes and R2 values are shown for linear regressions forced through the origin. or when jjAjjH + lEjj > 250 W m2, where H and lE are the measured sensible and latent heat fluxes. Values outside this broad range arise either because of instrumentation problems or because the site is subject to considerable advection. Wilson et al. [2002] found that H + lE < A at many eddy flux sites globally. This is also the case for most of the15 flux sites used in this study, as seen in the scatterplots in Figure 4 of daily average H + lE versus daily average A. The exceptions are Bondville, Niwot Ridge, Tumbarumba and Virginia Park, for which linear regression slopes are close to unity, with the worst cases being Hess and Kendall where the slopes are 0.70. We have assumed that A = H + lE to ensure internal consistency in the analysis used to estimate the model parameters for each site. Using the measured A and lE to estimate the parameters in the model for Gs would be incorrect since this assumes the lack of energy closure is due to errors in H alone. An alternative is to assume that the mean Bowen ratio (H/lE) is correct so that we can divide both H and lE by the slopes of the respective regression plots shown in Figure 4. A comparison of results for the two approaches are discussed later for the Hess and Kendall sites. [21] Aerodynamic conductance was calculated using Ga ¼

k 2 um ; ln½ðzm  d Þ=zom ln½ðzm  d Þ=zov

ð13Þ

where zm is the height of wind speed and humidity measurements, d zero plane displacement height, zom and zov are the roughness lengths governing transfer of momentum and water vapor, k von Karman’s constant (0.41), and uz is wind speed at height zm [Monteith and Unsworth, 1990]. The quantities d, z0m and z0v were estimated using d = 2h/3, z0m = 0.123h and z0v = 0.1z0m, where h is canopy height [Allen et al., 1998]. Atmospheric stability modifies the values of Ga calculated using equation (13) by up to ±25% but we have neglected stability effects because evaporation from dry canopies is relatively insensitive to errors in Ga. Similar arguments apply to uncertainties in the ratio z0v/z0m. The sensitivity of E to Ga is particularly weak at the daily time steps used in this study and this fact is exploited in the companion paper by Zhang et al. [2008], who assign constant values of Ga = 0.033, 0.0125 and 0.010 m/s for forests, shrubs, grassland and crops, respectively, because of the lack of routinely available, quality wind speed data at 1-km resolution. [22] Meteorological variables measured at each flux station have been used in this paper rather than large-scale meteorology that would likely be used in any operational application of the algorithms developed in this paper. This is not a significant weakness since both Cleugh et al. [2007] and Mu et al. [2007] showed that there is very little degradation in the performance of the PM equation when

7 of 17

W10419

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

W10419

Figure 5. Values of gsx and f extracted from Tables 2 and 3 in rank order for the various biomes for (a and b) the six-parameter model and (c and d) the two-parameter model (gsx and f). Biome codes are given in Table 1. meteorological data at 0.05° (latitude and longitude) scales are substituted for locally measured inputs. 3.2. Leaf Area Index [23] The Lai data needed to compute Gc, and hence Gs, were extracted for a 7 7 km2 area centered on each flux tower from the 8-day standard MOD15A2 collection 5 product [Myneni et al., 2002] from the ORNL-DAAC. These data were derived from MODIS, the Moderate Resolution Imaging Spectrometer mounted on the polarorbiting Terra satellite, which has a daily overpass at around 1030 local time. The MOD15A2 product contains four data layers: the fraction of photosynthetically absorbed radiation, fPAR and Lai, plus their respective quality control layers. The quality assessment (QA) flags in the database were used to check the quality of the MODIS-Lai data. At each grid, all poor quality Lai data were deleted and replaced by interpolated values obtained from a piecewise cubic, hermite interpolating polynomial [Zhang et al., 2006]. Then, all of the quality-controlled Lai data for each pixel were smoothed by the Savitzky-Golay filtering method that is widely used for filtering MODIS-Lai [Fang et al., 2008] and other remote sensing data.[Jonsson and Eklundh, 2004; Ruffin et al., 2008; Tsai and Philpot, 1998]. The MODIS Lai provides reasonable estimates for most Australian vegetation types except for the open forest and woodlands in eastern Australia [Hill et al., 2006]. [24] Figure 3b shows there is essentially no correlation between annual precipitation, P, and annual mean Lai for the 15 selected sites. Lai ranged from 3.6 for an evergreen broad

leaf forest with P = 1600 mm to Lai = 0.7 for a savanna with 400 mm of rainfall.

4. Results 4.1. Parameter Estimation [25] The generalized pattern search algorithm in MATLAB1 (The MathWorks, Inc.) was used to optimize for each site the parameters gsx, Q50, D50, kQ, kA and f in the model for Gs (equations (12) and (6)). This was done by minimizing the cost function F N  X 2 Emeas;j  ERS;j



j¼1 N  X

Emeas;j  E meas

2

;

ð14Þ

j¼1

where Emeas,j and ERS,j are the jth measured and simulated daily average evaporation rates, E meas is the arithmetic mean of the measurements and N is the number of samples. ERS,j was calculated by substituting the modeled Gs into equation (1), in conjunction with the required meteorological variables as measured at each flux station. [26] Optimized values for the six parameters in the surface conductance model are listed for the 15 sites in Table 2, with the last row listing the ranges allowed for each parameter in the optimization calculations. Figure 5a shows values of gsx extracted from Table 2 in rank order for the

8 of 17

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

W10419

W10419

Table 2. Site Name, Biome Code, and Optimized Values for All Six Parameters of the Surface Conductance Model for the 15 Flux Stations Used in the Analysisa Site Name

Code

gsx (m/s)

f

Q50 (W/m2)

D50 (kPa)

kQ

kA

a

Tonzi Hesse Virginia Park Howland Morgan SSA_OBS Kendall Hainich Bondville Mize Tumbarumba Santerem Km83 Griffin Niwot Ridge Mer Bleue Parameter ranges

wsa dec wsa dec dec con gra dec cer con ebf ebf con con wet –

0.0030 0.0067 0.0063 0.0038 0.0036 0.0020 0.0048 0.0046 0.0053 0.0039 0.0047 0.0076 0.0085 0.0027 0.0028 0.002 – 0.015

0.05 0.06 0.09 0.16 0.21 0.51 0.53 0.76 0.80 0.92 1.00 1.00 1.00 1.00 1.00 0.05 – 1.00

20.0 50.0 20.0 50.0 30.6 20.0 20.0 20.0 20.0 48.8 38.1 20.9 20.0 20.0 20.0 20 – 50

0.70 0.70 0.70 0.74 0.70 0.70 0.70 0.70 0.70 0.70 0.70 0.70 0.70 0.70 0.70 0.7 – 1.5

0.30 0.30 0.30 0.62 0.30 0.74 0.30 1.00 0.30 0.30 0.68 1.00 1.00 0.30 0.30 0.3 – 1.0

0.80 0.80 0.80 0.80 0.80 0.50 0.50 0.60 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.5 – 0.8

0.52 0.87 0.36 0.76 0.71 0.86 0.88 0.91 0.88 0.78 0.92 0.71 0.84 0.60 0.86 –

b (mm/d) 0.51 0.16 0.80 0.37 0.29 0.07 0.13 0.05 0.12 0.51 0.09 0.87 0.00 0.56 0.13 –

R2 0.52 0.87 0.48 0.81 0.73 0.82 0.84 0.86 0.81 0.80 0.84 0.74 0.69 0.43 0.79 0.73

Emeas Percent Error Percent Error (mm/d) Systematic Unsystematic 1.09 1.10 1.20 1.23 1.22 0.74 1.68 1.13 2.18 2.37 2.11 3.06 1.61 1.79 1.50 1.60

37.9 14.7 40.3 22.5 33.6 13.1 8.1 7.9 8.6 10.9 6.2 7.9 18.9 24.2 13.5 17.9

a

40.3 36.4 23.0 34.0 49.2 37.0 25.1 27.9 25.4 20.0 24.2 11.6 41.1 38.9 35.2 31.3

Also shown are the slope, a, and intercept, b, of the linear regression ERS6 = aEmeas + b, the R2 value, mean annual Emeas, and the percentage systematic and unsystematic root-mean-square error in average evaporation. Allowable parameter ranges used in the optimization are shown in the last row, as well as averages for Emeas and the percentage errors across all sites. The flux sites have been sorted according to increasing f value.

Figure 6. Time series of (top) 8-day averages for Emeas, ERS6, and Eeq and (middle) precipitation (cm/d) and Lai and (bottom) scatterplots of ERS6 versus Emeas for the Hesse, Tumbarumba, and Bondville sites. 9 of 17

W10419

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

W10419

Figure 7. Time series of 8-day averages for (top) Emeas, ERS6, and Eeq and (middle) precipitation (cm/d) and Lai and (bottom) scatterplots of ERS6 versus Emeas at the SSA Old Black Spruce site and two savanna sites, Virginia Park and Tonzi. various biomes. Maximum stomatal conductance ranged between 0.0020 and 0.0085 m/s, which is below the maximum of 0.012 m/s reported by K95 for various vegetation types, but is similar to the range reported by Isaac et al. [2004] for water stressed to well-watered crops and pastures. The soil evaporation parameter f reduces soil evaporation below the equilibrium rate when water availability limits evaporation. While this is a rather crude parameterization, the optimized values of f were reasonable, ranging from a low of 30% for the two savanna sites and for Morgan, 24.2% at Niwot Ridge, and an average of 13.8% for the other 11 sites. There is a significant inverse relationship between f and the systematic component of the error but not in the unsystematic component, 100eRMSE,u/Emeas, which has an average of 31.3%. [34] The scatterplot in Figure 8a shows daily average ERS,6 versus Emeas for all sites combined and linear regression yields a slope of 0.85, intercept of 0.21 mm/d and R2 = 0.82, N = 10630. The results show that good estimates of daily average evaporation are obtained when the PM equation is used with the model for Gs when all six parameters are optimized at each site and with MODIS Lai data. 4.3. Comparison With Cleugh et al. [2007] [35] Cleugh et al. [2007] obtained good agreement between measured evaporation rates at two Australian sites, Tumbarumba and at Virginia Park, when they used the simple linear model Gs = cLLai + Gs,min in the PM equation,

11 of 17

W10419

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

W10419

Figure 8. Predicted daily average daily evaporation rates (mm/d) plotted against average daily measured evaporation rates. Predictions use (a) the six-parameter model for Gs; (b) Gs = cL Lai, [Cleugh et al., 2007]; (c) equilibrium evaporation, Eeq; (d) the two-parameter model for Gs; and (e) the oneparameter model with gsx assigned from K95.

where cL is an empirical coefficient and Gs,min is the surface conductance controlling soil evaporation (set to zero in their analysis).The current data set was used to optimize cL for each site before calculating daily average evaporation rates, EcL, using the PM equation. Values of cL ranged from 0.0005 to 0.0037 (Table 4), compared to values of 0.0019 at Tumbarumba and 0.0025 at Virginia Park reported by Cleugh et al. [2007]. Figure 8b shows that daily average evaporation rates calculated using the simple linear model for Gs captures much of the variation in the measured evaporation, but the larger scatter (R2 = 0.61) indicates that the results obtained using EcL are less satisfactory than ERS6. Both models are clearly superior to using Eeq to estimate actual evaporation (Figure 8c), since Eeq is generally greater than the measurements, often by a factor of two or more, the

intercept of the regression is significantly greater than zero and the scatter in Eeq is unacceptably large (R2 = 0.57). 4.4. Reducing the Number of Free Parameters in the Gs Model [36] The general utility of the Gs model for estimating evaporation will be greatly enhanced if we can reduce the number of parameters to be estimated without substantially degrading model performance. To examine this possibility we explore two different approaches. The first is to take advantage of the observation, made earlier when discussing Figure 1, that Gs is relatively insensitive to Q50, D50 kQ and kA, and hence the optimization for gsx and f was repeated while holding the other parameters constant. The resultant two-parameter model for Gs is then used to calculate

12 of 17

LEUNING ET AL.: MODIS-LAI-BASED EVAPORATION MODEL

W10419

W10419

Table 3. Site Name, Biome Code, and Optimized Values of gsx and f of the Two-Parameter Surface Conductance Model for the 15 Flux Stations Used in the Analysisa Two-Parameter Optimization Site Name

Code

Tonzi Hesse Virginia Park Howland Morgan SSA_OBS Kendall Hainich Bondville Mize Tumbarumba Santerem Km83 Griffin Niwot Ridge Mer Bleue Average

wsa dec wsa dec dec con gra dec cer con ebf ebf con con wet

gsx (m/s) 0.0037 0.0063 0.0069 0.0029 0.0050 0.0022 0.0060 0.0040 0.0069 0.0044 0.0042 0.0066 0.0086 0.0038 0.0043 0.0051

f

a

b (mm/d)

0.05 0.13 0.09 0.25 0.05 0.55 0.50 0.84 0.79 0.70 1.00 1.00 1.00 1.00 1.00

0.48 0.84 0.34 0.73 0.67 0.86 0.88 0.91 0.87 0.75 0.90 0.74 0.82 0.59 0.85

0.54 0.22 0.82 0.42 0.33 0.07 0.12 0.04 0.12 0.59 0.13 0.73 0.02 0.55 0.15

R2

Emeas (mm/d)

Percent Error Systematic

Percent Error Unsystematic

0.47 0.87 0.45 0.81 0.69 0.81 0.83 0.86 0.79 0.77 0.83 0.65 0.67 0.41 0.76 0.71

1.09 1.10 1.20 1.23 1.22 0.74 1.68 1.13 2.18 2.37 2.11 3.06 1.61 1.79 1.50 1.60

41.0 17.8 41.3 25.5 38.1 13.6 8.4 8.1 9.4 12.6 7.1 7.4 20.7 24.9 14.4 19.3

40.9 36.0 23.1 32.4 50.5 37.9 26.2 28.5 26.6 20.5 24.5 14.8 41.7 39.6 37.6 32.1

a Also presented are the slope, a, and intercept, b, of the linear regression ERS2 = aEmeas + b, the R2 value, mean annual Emeas, and the percentage systematic and unsystematic root-mean-square error in average evaporation. Fixed parameter values are kQ = kA = 0.6, Q50 = 30 W/m2, and D50 = 0.7 kPa.

evaporation rates, defined as ERS2. The fixed values chosen for parameters Q50 = 30 W m2 and D50 = 0.7 kPa were guided by the results in Table 2. In principle, the extinction coefficient for visible radiation, kQ is expected to differ from that for available energy, kA, which accounts for the transfer of the sum of visible, near infrared and net thermal radiation. This is because reflection and transmission coefficients of leaves are considerably lower in the visible waveband than in the near infrared while in the thermal waveband the leaves have very high emissivity and low transmissivity. However, Gs is insensitive to values of the extinction coefficients, and a constant kQ = kA = 0.6 was used in the subsequent analysis. [37] Table 3 and Figures 5c and 5d show that optimizing just two parameters caused only minor changes in the values and relative ranking of gsx and f compared to those shown in Table 2 and Figures 5a and 5b for the six-parameter model. This provides further evidence that modeled evaporation rates are relatively insensitive to Q50, D50 kQ and kA.

Comparison of Tables 2 and 3 and Figures 8a and 8d shows that the two-parameter model for Gs performs almost as well as the six-parameter model in describing the variation in daily average evaporation rates across all sites. The linear regression of ERS2 versus Emeas has a slope of 0.83, intercept = 0.22 mm/d and R2 = 0.80. The nonzero intercept is largely due to the systematic bias of the model at the two savanna sites, Tonzi and Virginia Park. The model performs least well for the three sites with values of f 0.09, indicating the importance of soil evaporation at these sites and the desirability of finding a means to estimate the temporal and spatial variation of f using remote sensing rather than treating f as a fixed parameter. The average systematic root-mean-square error in daily mean E across all 15 sites was 0.31 mm/d (19.3% of 1.60 mm/d), with a range of 0.05 – 0.50 mm/d. Compared to the six-parameter model, the two-parameter model resulted in a