Two Improvements of an Operational Two-Layer Model for Terrestrial Surface Heat Flux Retrieval

Abstract

In order to make the prediction of land surface heat fluxes more robust, two improvements were made to an operational two-layer model proposed previously by Zhang. These improvements are: 1) a surface energy balance method is used to determine the theoretical boundary lines (namely ‘true wet/cool edge’ and ‘true dry/warm edge’ in the trapezoid) in the scatter plot for the surface temperature versus the fractional vegetation cover in mixed pixels; 2) a new assumption that the slope of the T_m – f curves is mainly controlled by soil water content is introduced. The variables required by the improved method include near surface vapor pressure, air temperature, surface resistance, aerodynamic resistance, fractional vegetation cover, surface temperature and net radiation. The model predictions from the improved model were assessed in this study by in situ measurements, which show that the total latent heat flux from the soil and vegetation are in close agreement with the in situ measurement with an RMSE (Root Mean Square Error) ranging from 30 w/m²∼50 w/m², which is consistent with the site scale measurement of latent heat flux. Because soil evaporation and vegetation transpiration are not measured separately from the field site, in situ measured CO₂ flux is used to examine the modeled λE_veg. Similar trends of seasonal variations of vegetation were found for the canopy transpiration retrievals and in situ CO₂ flux measurements. The above differences are mainly caused by 1) the scale disparity between the field measurement and the MODIS observation; 2) the non-closure problem of the surface energy balance from the surface fluxes observations themselves. The improved method was successfully used to predict the component surface heat fluxes from the soil and vegetation and it provides a promising approach to study the canopy transpiration and the soil evaporation quantitatively during the rapid growing season of winter wheat in northern China.

1. Introduction

Evapotranspiration (λE, soil evaporation and vegetation transpiration) from the land surface is an important link between the surface energy balance and the hydrologic cycle. Its accurate characterization is therefore very important in the study of the terrestrial ecosystem, climate dynamics and hydrologic cycle. At present, estimate of regional evapotranspiration has been made possible by using the remote sensing observations in combination with the surface meteorological data. In the past years, several remote sensing methods were developed to simulate surface-atmosphere interactions and to retrieve the terrestrial evapotranspiration over a wide range of spatial scales [1].

By treating the soil-vegetation system as a single uniform leaf, the big-leaf model simplified the mechanism of the energy exchange between the surface and the atmosphere, and therefore the regional scale evapotranspiration simulation is made. This category of models is simple and convenient to use, but the limitation is that this big-leaf approximation in the model is not applicable to surfaces with highly spatial heterogeneity due to large differences of surface energy exchange between soil and vegetation, such as in arid or semi-arid areas. Therefore, a two-layer model is proposed and the surface available energy is partitioned between soil and vegetation to overcome the limitation of the big-leaf model. These models have an improvement over the big-leaf models when applied to sparsely vegetated surfaces [2]. Another motivation to study the two-layer model is that China exhibits a highly heterogeneous land cover due to a large population and scattered residential areas. The two-layer model is more suitable in such areas.

In the existing two-layer models, the cores of the algorithm primarily lie in two aspects: (1) accurately decomposing surface temperature of mixed pixel (T_m) into soil temperature (T_soil) and vegetation temperature (T_veg); (2) obtaining accurate surface resistances, such as aerodynamic resistance, canopy resistance, residual resistance. In recent years, many attempts had been made to investigate the two issues. For instance, Norman and Kustas [3,4,5] used remote measurements of surface directional brightness temperature and some ancillary data to obtain soil temperature and vegetation temperature (called multi-angle method), applied Beer's law to partition net radiation of mixed pixel (R_n) and employed Monin-Obukhov similarity theory to compute aerodynamic resistance (r_a); Zhang et al. [6] presented a PCACA algorithm (PCACA, Pixel Component Arranging and Comparing Algorithm) and a layered energy-separating algorithm on the basis of triangle method and Bowen-ratio energy balance method to partition surface temperature, surface albedo (α_m) and net radiation of mixed pixel, and finally to estimate soil evaporation (λE_soil) and vegetation transpiration (λE_veg). Because the multi-angle satellite data is not always available, multi-angle method of surface temperature decomposing is limited for applications, In contrast, PCACA algorithm is more convenient because only single angle remote sensed data are required and it can be provided from most of the satellite data. Additionally, by using the layered energy-separating algorithm the core of which is Bowen-ratio energy balance method, the uncertainties in surface energy partitioning based on the Beer's law are reduced.

PCACA algorithm and layered energy-separating algorithm utilize the scatter plot of the surface temperature against vegetation fraction cover (T_m – f space) to determine the soil water status, like the trapezoid method. On the basis of: (1) the assumption that the configuration of T_m – f space is not primarily caused by differences in atmospheric conditions and soil attributes (eg. air temperature T_a, aerodynamic resistance r_a, surface reflectivity α) but by the variations of soil water availability; (2) the fact that iso-lines of equal soil water availability are nearly straight in T_m – f space, which was reported in previous studies on trapezoid method [7-11], Zhang et al. [6] indicated that T_soil values for all pixels at an iso-line are equivalent, so are for T_veg values, which is just like the case that while measuring the same area constituted by soil and vegetation at varying view angles, T_soil and T_veg are invariable and thus T_m observations only vary with vegetation fraction cover. Under this condition, T_soil and T_veg could be obtained by calculating the slopes of these iso-lines of equal soil water availability (dT_m/df). Detailed descriptions about the two algorithms and the trapezoid method can be found in Zhang et al. [6] and in Carlson et al. [9, 12], respectively. In this approach, an assumption of a uniform atmospheric environment and homogeneous soil surface is required. In most cases, however, this strict requirement can't be satisfied, especially on a regional scale. In addition, the identification of the trapezoid shape of T_m – f space in the trapezoid method requires a flat surface and a large number of pixels over an area with a wide range of soil wetness and vegetation fraction cover, which usually cannot be satisfied within a limited study area, and therefore some subjectivity would be introduced in the determination of the trapezoid shape and it will inevitably introduce some extra errors in theT_soil and T_veg calculations, The temperature separation finally influences λE estimate.

In this paper, to improve the accuracy of λE estimate using the two-layer model presented by Zhang et al. [6], two modifications were made to the PCACA algorithm and layered energy-separating algorithm mentioned above: (1) to ensure that the assumption that the configuration of T_m – f space is mainly controlled by soil water availability is reasonable, the effects of atmospheric conditions on surface temperature would be ignored by using the averaging method; (2) to identify the shape of the trapezoid bounded by true wet/cool edge' and ‘true dry/warm edge, a general method based on surface energy balance was used. Finally, comparisons between λE retrievals from the original model and the improved one, and in situ λE measurements were done to assess the improved method.

2. Model Descriptions

The two-layer model used in this study was presented by Zhang et al. [6], and in it the PCACA algorithm and layered energy-separating algorithm are the key algorithms. In the model, land cover is simplified as a mixture of two elements, namely, vegetation and bare soil. The energy fluxes are partitioned between the soil and vegetation, and energy exchange between vegetation and bare soil is negligible. Two parameters, albedo and surface temperature, are the main different characteristics of vegetation and bare soil, and lead to different interactions between them and atmosphere. Air temperature, air humidity, wind speed and aerodynamic resistance are approximated the same for vegetation and bare soil in the same pixel due to intensively atmospheric blending effect.

2.1 An interpretation of T_m – f space

Figure 1 provides a conceptual illustration of T_m – f space, where “true wet/cool edge” of the trapezoid is related to surface conditions of potential evapotranspiration and has minimum surface resistance to evapotranspiration (r_smin).

Soil water content equals to field capacity; ‘true dry/warm edge’ represents zero evapotranspiration and has maximum surface resistance to evapotranspiration (r_smax). If the positions of the two edges are determined, the shape and the structure of the trapezoid can be fixed and the consequent calculations of the surface heat fluxes could be done. To make the illustration easier to follow, four points are defined: T_sd and T_sw represent the points of true dry bare soil and true water saturated bare soil, respectively. T_vd represents true dry full-cover vegetation and T_vw represents true water saturated, full-cover vegetation. The above definitions all correspond to the ideal surface conditions, namely there exist driest and wettest bare soil and full-cover vegetation. In reality there are always insufficient number of pixels that can cover all kinds of soil wetness and vegetation fraction cover within the study area, which leads to a difficulty in determining “true wet/cool edge” and “true dry/warm edge”, as a result, “observed wet/cool edge” and “observed dry/warm edge” (dashed lines) are often defined according to the envelop shape of the actual scatter plot to represent actual two extreme soil moisture conditions and are used to replace “true wet/cool edge” and “true dry/warm edge” in most applications, although some errors would be introduced. Iso-line of equal vegetation fraction cover intersects “true dry edge” and “true wet edge” at true maximum temperature and true minimum temperature denoted as T_mi,max and T_mi,min, respectively. It has to be noted that for the trapezoid constructed by data with coarser pixel resolution (eg. 1km), surface temperature at “true dry edge” generally is higher than that at “observed dry edge”, contrarily surface temperature at “true wet edge” is lower than that at “observed wet edge” according to the findings of Carlson [9] about the scale issues of trapezoid space. The major reason is the absence of the driest and wettest pixels due to the great mixing effect of the pixel (vegetation/bare soil).

2.2. PCACA Algorithm

The PCACA algorithm is a method to decompose the mixed surface temperature. As mentioned above, its physical basis is the fact that soil water status can be represented by the configuration of T_m-f space, and like the triangle method, the crucial assumption is that surface temperature is mainly controlled by soil water availability [9, 10, 12].

The basic formulations used to decompose mixed surface temperature are shown in Eq. (1) and Eq. (2). The derivations of Eq. (2) can be found in Appendix I.

$$\sigma {\varepsilon }_m{T}_m^4=\sigma {\varepsilon }_{\mathit{\text{veg}}}f{T}_{\mathit{\text{veg}}}^4+\sigma {\varepsilon }_{\mathit{\text{soil}}}\left(1-f\right)T_{\mathit{\text{soil}}}^4$$

(1)

$$k=\frac{d{T}_m}{df}\approx T_{\mathit{\text{veg}}}-T_{\mathit{\text{soil}}}$$

(2)

where, ε_m, ε_veg and ε_soil are the broadband emissivities of mixed pixel, vegetation and bare soil, respectively; σ is the Stefan-Boltzmann constant; dT_m/df represents the slope of iso-line of equal soil water availability, expressed as k in the remainder of the paper. In the study, constant ε_veg of 0.97 and ε_soil of 0.95 were used. By simply weighting the fractional cover for vegetation and bare soil, mixed pixel emissivity ε_m for each pixel can be computed following Eq. (3) [13].

$${\varepsilon }_m=f{\varepsilon }_{\mathit{\text{veg}}}+\left(1-f\right){\varepsilon }_{\mathit{\text{soil}}}$$

(3)

From the above three equations, we can see that to solve T_soil and T_veg, k for each pixel is the key parameter and needs to be obtained firstly.

According to Figure 1, two procedures were needed to derive k in the study: (1) determining the shape of the trapezoid bounded by ‘true dry edge’ and ‘true wet edge’, namely locating the ‘true dry edge’ and ‘true wet edge’ on the trapezoid, and thereby calculating the two slopes of the two bounder lines (k_u is for upper bound, k_L is for lower bound). A detailed method for this will be illustrated in section 3; (2) linearly interpolating k between the highest and the lowest surface temperature to obtain the slope for each pixel based on the equation below:

$$k_i=\frac{T_{mi}-T_{mi,min}}{T_{mi,max}-T_{mi,min}}\left(k_u-k_L\right)+k_L$$

(4)

Considering that the precise formulation of the relationship is actually unknown, linear interpolation is a reasonable approximation according to the study on the configuration of T_m – f space by Calson [9, 10] and Moran [14].

2.3 Layered Energy-separating Algorithm

Essentially, the aim of the Layered Energy-separating algorithm is to calculate Bowen-ratio (β, the radio of sensible heat flux to latent heat flux) of soil and vegetation expressed as β_soil and β_veg, respectively. By using the relationship between Water Deficit Index (WDI), evapotranspiration and potential evapotranspiration (λE₀) illustrated by Moran et al. [14], Eq. (5) can be derived:

$${\beta }_i\approx \frac{T_{m,max}-T_{m,min}}{T_{m,max}-T_{mi}}-1$$

(5)

where subscript i represents pixel i. From Eq.(5) we can see that β can be obtained conveniently from the data shown on Figure.1.

2.4 Estimation of Other Core Variables

After T_s, T_v, β_s and β_v are obtained using the above methods, net radiation at the soil surface and at the vegetation surface (R_n,soil, R_n,veg) can be calculated following Eqs. (6) and (7), and thereby soil heat flux and net radiation of mixed pixel can be estimated using the empirical formulation of Eq. (8) and Eq. (9), as used widely in previous studies [15-17].

$$R_{n,\mathit{\text{soil}}}=S_0\left(1-{\alpha }_{\mathit{\text{soil}}}\right)+\sigma {\varepsilon }_{\mathit{\text{sky}}}{T}_{\mathit{\text{sky}}}^4-\sigma {\varepsilon }_{\mathit{\text{soil}}}{T}_{\mathit{\text{soil}}}^4$$

(6)

$$R_{n,\mathit{\text{veg}}}=S_0\left(1-{\alpha }_{\mathit{\text{veg}}}\right)+\sigma {\varepsilon }_{\mathit{\text{sky}}}{T}_{\mathit{\text{sky}}}^4-\sigma {\varepsilon }_{\mathit{\text{veg}}}{T}_{\mathit{\text{veg}}}^4$$

(7)

$$G\approx 0.3\left(1-0.9f\right)R_{n,\mathit{\text{soil}}}$$

(8)

$$R_n=f{R}_{n,\mathit{\text{soil}}}+\left(1-f\right)R_{n,\mathit{\text{veg}}}$$

(9)

where S₀ is the solar incident total radiation and is regarded as spatially uniform for clear sky conditions at the regional scale, usually obtained from standard meteorological station; α_soil and α_veg are the albedo of bare soil and vegetation, also calculated by the PCACA algorithm (seen from Appendix И); ε_sky is the average sky emissivity and is approximately set to 1.0 in the study. T_sky is average sky temperature and usually approximates to the temperature at 37° view sky angle [18], detailed calculation about which will be described in the Section 3.

In terms of the Bowen-ratio energy balance method, soil evaporation (λE_soil) and vegetation transpiration (λE_veg) can be retrieved based on Eq. (10). When energy exchange between vegetation and bare soil is neglected, λE of a mixed pixel can be described as a linear combination of λE_soil and λE_veg, expressed as Eq. (11).

$$\lambda E_{\mathit{\text{soil}}}=\frac{R_{n,\mathit{\text{soil}}}-G}{1+{\beta }_{\mathit{\text{soil}}}},\lambda E_{\mathit{\text{veg}}}=\frac{R_{n,\mathit{\text{veg}}}}{1+{\beta }_{\mathit{\text{veg}}}}$$

(10)

$$\lambda E=f\lambda E_{\mathit{\text{veg}}}+(1-f)\lambda E_{\mathit{\text{soil}}}$$

(11)

3. Two improvements for the two-layer model

3.1 locating the true dry edge in T_m – f space

As mentioned above, the locations of true dry edge and the true wet edge are crucial in the application of PCACA algorithm. The previously used method to determine the trapezoid boundary often leads to uncertainties because of subjective judgement. To reduce the errors from this respect, a physically based method, which takes account of the surface energy balance, is presented in this study.

According to Figure.1, the four corner points, T_sd, T_sw, T_vd and T_vw can determine the envelop shape of the trapezoid, that is to say, as long as their values are obtained. Consequently, the true dry edge and the true wet edge can be determined. In the study, surface energy balance method was adopted to compute their values. For pixels at the true dry edge, Eq.(12) is used because λE=0. Substituting R_n, G, H for Eq.(9), Eq.(8) and Eq.(13), we obtained Eq.(14) for T_sd. In the same way, T_vd is formulated as Eq.(15).

$$R_n-G=H$$

(12)

$$H=\frac{\rho Cp\left(T_s-T_a\right)}{r_a}$$

(13)

$$T_{sd}=\frac{0.7[S_0(1-{\alpha }_{sd})+\sigma {\varepsilon }_{\mathit{\text{sky}}}{T}_{\mathit{\text{sky}}}^4]+\frac{\rho Cp}{r_{\mathit{\text{sda}}}}{T}_{\mathit{\text{sda}}}}{\frac{\rho Cp}{r_{\mathit{\text{sda}}}}+0.7\sigma {\varepsilon }_{\mathit{\text{soil}}}{T}_{sd}^3}$$

(14)

$$T_{vd}=\frac{0.7[S_0(1-{\alpha }_{vd})+\sigma {\varepsilon }_{\mathit{\text{sky}}}{T}_{\mathit{\text{sky}}}^4]+\frac{\rho Cp}{r_{\mathit{\text{vda}}}}{T}_{\mathit{\text{vda}}}}{\frac{\rho Cp}{r_{\mathit{\text{vda}}}}+0.7\sigma {\varepsilon }_{\mathit{\text{veg}}}{T}_{vd}^3}$$

(15)

where the subscripts sd, vd represent true driest bare soil and true driest full-cover vegetation, respectively. ρ is density of air, c_p is the volumetric heat capacity of air, r_a is the air dynamic resistance, σ is the S-B coefficient, α_sd is the albedo of dry bare soil.

From these two equations, we can see that T_sd or T_vd could not be iteratively computed until parameters of α_sd, α_vd, T_sky, r_sda, r_vda, T_sda, T_vda are acquired. From the T_m – f space, the observed driest bare soil and the observed driest full-cover vegetation can be found. Although there are surface temperature differences between them and the true driest bare soil and true driest full-cover vegetation, they represent relative driest and wettest soil conditions. It means that aerodynamic resistance and air temperature of the observed driest bare soil and the observed driest full-cover vegetation are approximate to that of T_sd, T_vd points under the conditions that the variations in atmospheric conditions are very small due to the spatially uniformity. In the study, we chose 50 pixels around the upper-left corner in the trapezoid representing the observed driest bare soil and 50 pixels around upper-right corner in the trapezoid representing the observed driest full-cover vegetation, and the highest T_a and r_a in each 50 pixels were selected as T_sda and r_sda, T_vda and r_vda, respectively. The method of retrieving the spatial distribution of T_a and r_a will be described in the following. The calculation of r_a requires e_a and r_s, the retrievals about which will also be described. As for α_sd and α_vd, we selected the lowest albedo values of bare soil and full-cover vegetation within the whole scene of remotely sensed image (pixels) because low albedo would result in high surface temperature at the same soil water content, judging from Eq. (14) and Eq. (15). The item of $\sigma {\varepsilon }_{\mathit{\text{sky}}}{T}_{\mathit{\text{sky}}}^4$ represents downward long wave radiation that usually has small spatial variability for clear sky, thus can be obtained from the measurements of meteorological stations at the satellite overpass time. The following are the calculations of spatial distributions of T_a, e_a, r_s and r_a:

a) Estimation of air temperature (T_a)

Surface temperature, as a heat or cold source, influences the variations of air temperature by heating or cooling near-surface atmosphere. In most cases, high surface temperature is accompanied by high air temperature and low surface temperature is accompanied by low air temperature. By using this relationship between them and assuming the following ratio expressed in Eq. (16), an interpolation method is presented to map the distribution of air temperature. The numerator and the denominator both represent radiant emission. Because the effects of neighborhood pixels diminish as the distance to the pixel increases, a weight (w) derived from an inverse distance weight (IDW) method is assigned to each neighboring observations, and is given in Eq. (17).

$$\frac{\sum _{i=1}^n{w}_i\cdot \sigma \cdot {\varepsilon }_{ai}\cdot T_{ai}^4}{\sum _{i=1}^n{w}_i\cdot \sigma \cdot {\varepsilon }_{mi}\cdot T_{mi}^4}=\frac{\sigma \cdot {\varepsilon }_{aj}\cdot T_{aj}^4}{\sigma \cdot {\varepsilon }_{mj}\cdot T_{mj}^4}$$

(16)

$$\begin{array}{cc}{w}_i=\frac{d_i^{-p}}{\sum _{i=1}^n{d}_i^{-p}},& {\varepsilon }_a=1.24\cdot {(\frac{e_a}{T_a})}^{1/7}\end{array}$$

(17)

where i represents the pixel where air temperature measurement is made, j indicates the pixel where no measurements are available and estimate is required; n is the number of air temperature observations; d is the distance between pixel i and pixel j; p is an exponent. The higher the exponent is, the larger the influence of the closest observations on estimate is; p is set to 2 in the study; σ is the Stephen Boltzmann constant; ε_a is the emissivity of air calculated from Eq. (17) [19].

b) Estimation of actual vapor pressure near surface (e_a)

Like the relationship between air temperature and surface temperature, there are also strong interactions between near-surface actual vapor pressure and soil moisture. If there is no horizontal and vertical advection, vapor in the atmosphere mainly comes from soil water by an evapotranspiration process. Contrarily, the vapor gradient between surface and air influences the intensity of evapotranspiration. By assuming the following ratio relationship between near-surface actual vapor pressure and soil water, we interpolated e_a, where soil moisture status was characterized quantitatively by the soil apparent thermal inertia:

$$\frac{Rn}{2\left(T_s-T_{min}\right)\Delta t^{\frac{1}{2}}}\cdot \frac{\sum _{i=1}^n{w}_i\cdot e_{ai}}{\sum _{i=1}^n{w}_i\frac{R_{ni}}{\left(T_{mi}-T_{min}\right)}}=\frac{e_{aj}}{\frac{R_{nj}}{\left(T_{mj}-T_{min}\right)}}$$

(18)

where i, j, w_i, T_s have the same meanings as in Eqs. (16) and (17); T_min is the minimum surface temperature during the daytime which usually occurs before sunrise when R_n=0 and for all pixels it can be assumed to be the same value [20]. Δt is the time interval between sunrise and the satellite overpass time.

c) Estimation of surface resistance to evapotranspiration (r_s)

In the retrieval of evapotranspiration, r_s is used to correct the difference between the vapor pressure at the surface (e_s) and the saturated vapor pressure at the evaporating front (e_s*). In theory, r_s ranges from ∞ to 0 corresponding to the surface conditions of potential evapotranspiration and zero evapotranspiration, respectively, namely the true dry edge and the true wet edge. However, in reality the condition of zero evapotranspiration (r_smax = ∞) rarely occurs for vegetated surface even in semi-arid environment primary due to root zone soil water uptake, consequently, we selected a pixel closest to the observed driest bare soil where a meteorological station is located to calculate r_smax by λE, e_a, u measurements at the satellite overpass time. r_smax is about 1000 (s m^-1) according to the calculation, which is in agreement with Qiu's observations [21]; r_smin is set to 0 in the study.

After the upper and lower limits of r_s are determined, r_s is interpolated linearly within each f interval between the lowest and the highest temperature. Following the illustrations in Figure 1, r_si for a pixel at (f_i , T_i) equals:

$$r_{si}=\frac{T_{si}-T_{si,min}}{T_{si,max}-T_{si,min}}\left(r_{smax}-r_{smin}\right)+r_{smin}$$

(19)

The interpretation is similar to the method used by Stisen et al. [11]. The difference is that Stisen et al. used Priestly-Taylor parameter to represent an effective surface resistance to evapotranspiration.

d) Estimation of aerodynamic resistance (r_a)

Besides air temperature, aerodynamic resistance is a site-specific variable and can not be retrieved directly by remote sensing. Although Monin-Obukhov similarity theory has been widely used to estimate it, the accurate calculations for the spatial distributions of roughness length (Z₀), wind speed (u) and the atmospheric stability parameters are very difficult. In this study, we adopted the energy balance method.

An expression of r_a is obtained on the basis of the energy-balance by substituting Rn, G, H and λE for Eq. (9), Eq. (8), Eq. (13) and Eq. (20):

$$\lambda E=\frac{\rho c_p\left(e_s^{*}-e_a\right)}{\gamma \left(r_s+r_a\right)}$$

(20)

$$\begin{array}{l}{a}_1{r}_a^2+a_2{r}_a+a_3=0\hfill \\ a_1=(0.7+0.27f)\cdot \gamma \cdot R_n\hfill \\ a_2=(0.7+0.27f)\gamma \cdot R_n\cdot r_s-\rho c_p\gamma (T_m-T_a)-\rho c_p(e_s^{*}-e_a)\hfill \\ a_3=-\rho c_p\gamma (T_m-T_a)\cdot r_s\hfill \end{array}$$

(21)

where γ is the psychometric constant, e_s* is the saturated vapor pressure at T_m estimated using the classic formulation with regard to surface temperature, as Eq. (22), other parameters have the same meanings as the above illustrations. Here, T_a, e_a and r_s are estimated by the approaches mentioned previously.

$$e_s^{*}=6.1078exp(\frac{17.27\cdot T_m}{T_m+237.3})$$

(22)

Above all, all necessary parameters to calculate T_sd, T_vd can be retrieved according to the Eqs. (12) – (22) in combination with limited ancillary data. Here, T_m and α_m is obtained from MODIS (Moderate Resolution Imaging Spectroradiometer) standard land data products, MOD11 and MOD02 [22, 23], through the NASA Earth Observing System Data Gateway.

3.2 Locating the true wet edge in T_m – f space

Many studies [9, 24, 25] indicated that the dry edge is more evident than the wet edge in the T_m– f space. There often exist outlying points exhibiting a tail toward low values of temperature and f. Such points may represent anomalous surfaces mainly due to cloud contamination and usually are discarded from the analysis [10]. Thus, although T_sw and T_vw also can be parameterized on the basis of energy balance equation like T_sd and T_vd, the inputs (α_sd, α_vd, r_sda, r_vda, T_sda, T_vda) cannot be obtained by the above methods. Taking account of the fact that the surface radiant temperature of dense vegetation is very close to the ambient air temperature [9, 26, 27], we took the average air temperature at f =1 as T_vw. As for T_sw, we adopted an approximation of using the surface temperature of standing water body (such lake) as the surface temperature of true wet bare soil T_sw, that is to say, standing water body is regarded as the surface of potential evapotranspiration. In fact, it is not uncommon to find some patches of standing water body in a remote sensed image. In this study, the surface temperature of Dongping Lake (35.965°, 116.81°; water area: 209 km²) was used as T_sw. In applications, we found that the pixels with mixed surface temperatures below the true wet edge all scattered around the cloud pixels and the coastline pixels.

3.3 Physical illustrations for the uncertainties using the above locating methods

Using the above methods, the positions of the true dry edge and the true wet edge in the triangle can be located. When surface radiant temperature is mainly dominated by surface soil water content, the following relationship between surface temperature and soil water content is tenable (Figure 2), which has been illustrated in the Qiu's experiment [21]. Figure 2 suggests that when the soil water content ranges from 80% to 100% of the soil saturation, evapotranspiration will happen approximately at the level of potential Evapotranspiration. On the contrary, when the range of the soil water content is from the wilting point to the driest, evapotranspiration is close to zero, and in the two extreme conditions, surface temperature can be regarded as a constant due to the same evaporative cooling effect. On the basis of this, although the true driest and the true wettest points can't be determined using the above methods, the errors will be acceptable in the determination of T_sd, T_sw, T_vd and T_vw due to the constant surface temperature.

3.4 Elimination of the effects of atmospheric conditions on surface evapotranspiration

Same as the triangle method, an important assumption for the PCACA algorithm is that the surface evapotranspiration is primarily constrained by soil water availability, based on which T_soil values for all pixels with an equal water availability are identical, so are for T_veg values, that is to say, the configuration of T_s and f is mainly caused by the variation of the soil water availability and is irrelevant to the differences in atmospheric conditions and surface properties [11]. However, in reality it is not true that the atmospheric conditions and surface properties are entirely homogeneous within the large study area. There exist some differences eg. in wind speed, air temperature, surface roughness length and surface albedo, which are all site-specific variables. As a result, it will lead to the the variability of the surface evapotranspiration and thereby influence the surface temperature and the configuration of T_s and f. To satisfy the assumption, the effects of these factors need to be eliminated.

In the study, the effects of four controlling factors (r_a, e_a, T_a and albedo) on surface temperature were eliminated. According to the assumption, r_a, e_a, T_a and albedo values of each pixel should be equal, namely they are spatially homogeneous, to ensure that only water availability controls surface temperature. To meet this requirement, the average values of r_a, e_a, T_a and albedo in the image are assigned to each pixel and a new mixed surface temperature (T_mi’ ) is calculated based on energy balance equation. T_mi’ is the assumed temperature controlled only by soil water availability and therefore the new trapezoid constructed by T_m and f is more meaningful for the PCACA algorithm. It has to be noted that although the configuration of T_m’ – f space is different from that of the T_m – f space, the locations of the true dry edge and the true wet edge in the trapezoid don't change since the two theoretical extreme soil water status at the satellite overpass time are the same. On the basis of T_m’, T_veg’ (vegetation temperature controlled only by soil water availability) and T_soil’ (soil temperature controlled only by soil water availability) can be obtained. Now the problem is that the calculated T_veg’ and T_soil’ cannot represent the actual temperature of vegetation and soil after the above transformation., Therefore they can't be directly used in the consequent calculations. To estimate the true vegetation temperature and true soil temperature (T_veg and T_soil) from T_veg’ and T_soil’, the following method is used.

Assuming that the thermal energy fraction assigned to the vegetation and the soil are constants, the following expression [Eq. (23)] can be approximated. The left item represents the proportion of thermal energy difference for vegetation or soil to the total thermal energy difference, and the right item represents the proportion of vegetative or soil thermal energy to the total thermal energy. It suggests that changes of environmental conditions bring same effects on vegetation temperature, soil temperature and mixed surface temperature, for example, T_soil, T_veg and T_m values would all increase with the increase of solar radiation intensity and all influenced by wind speed. It is not likely to happen that the T_soil increases, but T_veg decreases in the same atmospheric conditions. Therefore, the assumption is reasonable in practice and it is different from Beer's law using fractional vegetation cover as the weight to partition thermal energy. Using this equation, T_soil⁴ - T_veg⁴ can be solved as Eq. (24).

$$\begin{array}{cc}\frac{T_{\mathit{\text{veg}}}^4-T_{\mathit{\text{veg}}}^{\prime 4}}{T_m^4-T_m^{\prime 4}}\approx \frac{T_{\mathit{\text{veg}}}^4}{T_m^4},& \frac{T_{\mathit{\text{soil}}}^4-T_{\mathit{\text{soil}}}^{\prime 4}}{T_m^4-T_m^{\prime 4}}\approx \frac{T_{\mathit{\text{soil}}}^4}{T_m^4}\end{array}$$

(23)

$$T_{\mathit{\text{veg}}}^4-T_{\mathit{\text{soil}}}^4=\frac{T_m^4}{T_m^{\prime 4}}(T_{\mathit{\text{veg}}}^{\prime 4}-T_{\mathit{\text{soil}}}^{\prime 4})$$

(24)

Combined with Eqs. (1) and (24), T_soil and T_veg can be solved. T_soil and T_veg can be directly obtained from Eq. (23), at the same time the uncertainties of the solution can be reduced to the minimum by this method because the difference between T_soil and T_veg is used instead of the absolute temperatures and Eq.(1) helps to maintain the thermal energy balance in the calculation of T_soil and T_veg values.

4. Study area and field measurements

The study area is located in the North China Plain and ranges from 35.2N to 40.84N in latitude, from 113.68E to 119.54E in longitude. The land use in the area is dominated by the rotating cropping of winter wheat and summer maize. Millet, soybean and cotton are also scattered planted in summer [28]. According to the traditional tillage practice, winter wheat is sown in early October, harvested in early or mid June next year, and summer maize is planted in early to mid June and harvested at the end of September. The soil is mostly silt, light loam and medium loam. Annual precipitation is about 600 mm, more than 50% of which falls during the summer monsoon between July and September. The groundwater table varies from 1.5 m to 3.5 m with an average of 2.5 m.

In the study, field measurements from 135 standard meteorological stations were used. The measurements include air temperature and actual vapor pressure at 2 m height above the surface, solar incoming radiation, surface radiative temperature, wind speed, upward longwave radiation, upward shortwave radiation, downward longwave radiation and downward shortwave radiation. Figure.3 shows the spatial distribution of these stations in the study area, where the triangle symbol (▲) marks the location of Yucheng Agro-ecosystem Station that can represents the largest agricultural area in the North China Plain [29]. Sensible heat flux and latent heat flux have been continuously measured by the Eddy Correlation (EC) system which is composed of a 3D sonic anemometer and an open path CO₂/H₂O analyzer since 2002. H, λE and CO₂ fluxes were originally sampled at 10 Hz and values averaged over 30 min were used in the study to validate the above mentioned methods. The rectangle symbol (■) shows the location of Dongping Lake. Water surface temperature has been measured for seven years since 2001 and was used to validate the MODIS land surface temperature products used in this study, and thereby to ensure the accuracy of T_sw values.

5. Satellite data

MODIS land data products, including MOD11 (Land Surface Temperature), MOD03 (Geolocation Data Set), MOD05 (Total Precipitable Water), MOD02 (Calibrated Geolocated Radiance) and MOD35 (Cloud Mask), for clear sky during springtime between Mar and June when winter wheat is the dominant crop were used in the study. The visible channels of MOD02 were processed with atmospheric correction using the SMAC algorithm [30] with the combination of the MOD03 and MOD05 data. The cloud was masked out based on MOD35 data. The basic variables used as inputs in the model consist of albedo, fractional vegetation cover and surface temperature, based on which the calculation of R_n, PCACA algorithm and Layered Energy-separating Algorithm were applied.

The instantaneous surface albedo was obtained by averaging reflectance values for several visible and near-infrared channels with the wavelength as the weight. This would introduce some errors because channel reflectance values observed by a satellite are reflectance at only one Sun-target-sensor configurations and it is generally different from the hemispherical reflectance. But we assumed its influence to be small as pointed by Nishida [27]. The fractional vegetation cover was estimated from the normalized difference vegetation index (NDVI) [27], as Eq.(25).

$$f=\frac{\mathit{\text{NDVI}}-{\mathit{\text{NDVI}}}_{min}}{{\mathit{\text{NDVI}}}_{max}-{\mathit{\text{NDVI}}}_{min}}$$

(25)

where NDVI_max and NDVI_min are NDVI values for full cover vegetation (f=1) and bare soil (f=0). According to the field NDVI measurements performed at Yucheng station, NDVI_min=0.09 and NDVI_max=0.78. Other variables, such as R_n, G, T_a, e_a, r_s, were obtained using the above mentioned methods.

6. Results and Discussion

Figure 4 shows the comparisons between T_soil, T_veg obtained from the original model and T_soil’, T_veg’ obtained from the improved model at an iso-line of equal water availability.

Apparently, the results of T_soil’, T_veg’ are in closer agreement with the above assumptions that soil surface temperature for all pixels at an iso-line are identical, so are the vegetative surface temperature. It provides a better physical foundation for the model.

Due to the absence separate measurement of the soil evaporation, soil heat flux, vegetation transpiration and vegetation heat flux from the field sit, estimates of R_n-G and λE from the model were compared to the measurements from Yucheng station, as shown in Figure 5 for available energy (R_n-G), in Figure 6 for λE.

Because direct comparison of λE_soil and λE_veg cannot be performed, the relationship between modeled λE_veg and measured CO₂ were used to indirectly validate the two-layer model since CO₂ fluxes are closely related to λE_veg[31, 32]. Figure 7 shows the scatter plot between modeled λE_veg and measured CO₂ flux. Note that the minus sign (-) for CO₂ fluxes means that the flux transferring direction is from up to down. The performance of the model was evaluated using the root mean squares difference (RMSD) and the mean absolute difference (MAD), which are defined as Eq.(26) and Eq.(27), respectively.

$$\mathit{\text{RMSD}}={[\sum _{i=1}^n{(P_i-O_i)}^2/n]}^{1/2}$$

(26)

$$\mathit{\text{MAD}}=\sum _{i=1}^n\left|P_i-O_i\right|/n]$$

(27)

where n is the number of observations, P represents the model estimated value, O represents the observed value.

The figures demonstrate that (1) no significant bias is found between the modeled and measured R_n-G with RMSD=46.3 and MAD=40.2. The high correlation coefficient (r²=0.79) suggests that R_{n, soil} and R_{n, veg} calculated from T_soil and T_veg are reasonable because Rn is calculated by them; (2) estimates of λE tends to overestimate for lower λE values and shows larger bias than R_n-G resulting in 54.1 of RMSD and 47.7 of MAD, which implies a larger uncertainty of the model in low λE conditions; (3) in general, seasonal variations in modeled λE_veg and measured CO₂ flux shows good agreement though few points in 2006 showed different variations perhaps influenced by horizontal/vertical advection or other factors. Vegetation transpiration generally increases with the crop growth, at the same time, since vegetation absorbs more CO₂ due to the more active photosynthesis in the daytime, CO₂ fluxes also increases. On the contrary, λE_veg and CO₂ fluxes both became very small after crop harvest or in winter. In the study, winter wheat was harvest in June, therefore λE_veg and CO₂ fluxes rapidly decreased on June 15 and June 19 in 2005. Figure 7 indirectly proves the soundness of the two-layer model to estimate the soil evaporation and vegetation transpiration separately.

In fact, there are two other factors which contribute to the uncertainty in the above validation. One is that point measurements usually can not represent the whole MODIS pixel (1 km²) because of the large scale disparity between them. The other is that the observed R_n, λE, H and G values at field cannot meet the surface energy balance closure resulting from many possibilities, such as instrumental errors, horizontal/vertical advections [33-35], however, the basis of the model used in the study is the surface energy balance, as most remote sensing models of evapotranspiration, such as SEBS [36], SEBAL [37], N95 [3]. Therefore, it is not uncommon to see the differences between modeled fluxes and observed ones. Compared with other studies, the RMSD of 54.1 (w/m²) for modeledλE are in a acceptable range [17, 27, 38, 39].

Figure 8 shows the λE_veg maps on May 2, 2005 for the North China Plain retrieved by the improved model and the original model, respectively. Differences were found between them. The λE_veg value from the improved model varies from 222 (w/m²) to 456 (w/m²), while the original λE_veg values varied from 50 (w/m²) to 560 (w/m²). In terms of the experiences and the knowledge on the vegetation transpiration, a small λE_veg dynamical range is more reasonable in the rapid growing season of winter wheat, that is to say, the results from the improved model is better than that of the original model although it is difficult to quantitatively evaluate the two results due to the absence of field measurements of the vegetation transpiration and soil evaporation at the satellite pixel scale. Furthermore, according to the above illustrations, the physical basis of the improved model is strengthened against the original model.

7. Conclusions

This paper presented two improvements of a two-layer model for the estimation of the land surface heat fluxes. The weakness in the original model was identified: (1) a subjective method to determine the true boundary lines from the scatter plot for the surface temperature of mixed pixel versus the fractional vegetation cover; (2) the assumption that the configuration of T_m – f space is mainly controlled by soil water availability, are not physically realistic. To investigate the two issues and obtain more realistic separation of the surface temperature for the soil and vegetation components, and thereby more accurate partitioning of the surface fluxes for the soil and vegetation components, introducing the surface energy balance method solves the problem (1) and the method of eliminating the effects of four controlling factors (r_a, e_a, T_a and albedo) on surface temperature solves the problem (2). At the same time, the interpolation methods of T_a and e_a, the methods for acquiring spatial distributions of r_a and r_s were also investigated. Finally, the improved model was applied to the North China Plain. The results showed good agreement with in situ terrestrial surface fluxes measurements. Furthermore, by comparing the seasonal variations of vegetation transpiration and CO₂ flux, and the vegetation transpirations retrieved by the improved model and the original model, the effectiveness of the improved two-layer model for estimating soil evaporation and vegetation transpiration were indirectly proved.