methodology

The flowchart shows the whole process from input to the creation of the final early warning products. It starts with the processing of images from METEOSAT geostationary satellite and WMO-GTS precipitation data by the energy balance monitoring and rainfall mapping systems. These two systems generate the energy and water balance products, which are then used as input for the desertification monitoring and crop growth simulation models. The final output, desertification indices and crop yield forecasts, mainly find their use in drought, food and famine early warning applications. Products generated by the energy balance monitoring and rainfall mapping systems are especially useful for hydrological and meteorological applications.

1) Input data

- Hourly received visible(VIS) and thermal infrared(TIR) images from the METEOSAT geostationary satellite. The visible band covers the wavelengths from 0.3 to 1.05 μm. The thermal infrared band covers the wavelengths from 10.0 to 13.1 μm.

- Precipitation data from the Global Telecom System of the World Meteorological Organisation (WMO-GTS)

2) Pre-processing

- Cloud level estimation using temperature thresholding

- Counting cloud level durations (CD)

- Extracting local noon and local midnight sub-frames

- Detection, statistics and repair of defect scanlines

3) ENERGY BALANCE MAPPING

Calibration

TIR channel : based on calibration coefficients provided with the satellite data

VIS channel : vicarious calibration using Cumulonimbus cloud albedo as reference

Calibration converts satellite VIS counts into planetary albedo (A') and satellite TIR counts into planetary temperature (T₀')

VIS atmospheric correction

Based on a two-flux global radiation transmission model as proposed by Kondratyev (1969). The model was extended to include both absorption and scattering. The effect of the atmosphere is parameterized in terms of optical depth (t or tau). Once the optical depth is known, the model allows to derive the atmospheric transmissivity (t) and the surface albedo (A) from the planetary albedo (A'). The fraction of solar radiation absorbed at the earth surface then follows from t*(1-A). The following figures show the surface albedo (A) and the absorbed radiation t*(1-A) as a function of the planetary albedo (A').

Atmospheric correction comes down to determining the optical depth (t). This is done using the darkest land pixels in the image as a reference. These represent dense forest with a typical surface albedo of 7%.

TIR atmospheric correction

The relation between the planetary temperature (T₀') and surface temperature (T₀) may be represented by:

k = cos i_m * (T₀-T_a) / (T₀'-T_a)

Here, k is an atmospheric correction coefficient and i_m is the satellite zenith angle. The air temperature T_a is also derived from the satellite data, as explained in the next section. The atmospheric correction comes down to finding a corresponding pair of T₀ and T₀' values and then calculating k. This is done by finding the driest pixels (highest temperatures), for which we assume LE= 0. In this specific case the actual surface temperature may be calculated from the net radiation (In) with T0=Ta+In/a, where a is the heat transfer coefficient.

Air temperature mapping

A plot of noon planetary temperatures (T₀'_n) versus midnight temperatures (T₀'_m) shows a linear relation. This relation is determined by means of regression analysis: T₀'_n= a.T₀'_m + b. Along this line heat exchange wit the atmosphere is variable. In case of perfect heat exchange there is no temperature difference between the surface and the atmosphere and consequently : T₀'_n = T₀'_m = T_a. Combining these two relations, the air temperature is found from the regression constants with T_a=b/(1-a). This air temperature is to be interpreted as the temperature at the top of the atmospheric boundary layer. The procedure is applied to a sub-window and the air temperature is assigned to the center pixel. Subsequently the window is shifted over the image and the procedure is repeated for each pixel. In this way a complete air temperature map is obtained.

Cloud detection

A cloud detection algorithm is used which discriminates clear and cloudy pixels. A minimum albedo map (Amin) is extracted from a sequence of 10 days noon visible images. A maximum noon surface temperature map (Tmax) is obtained by searching the highest temperature within a shifting window of 200*200 pixels. Both are considered to represent cloud free conditions. Subsequently several tests are carried out to determine if a pixel is cloudy or not. The most important ones are:

1) A' ³ A_min + DA

2) T₀' £ T_max - DT

DA and DT are thresholds, which have been determined empirically.

Global radiation

The transmissivity of the atmosphere (t) is derived together with the surface albedo as discussed in section 4. The global radiation at noon (I_gn) is then obtained with:

I_gn = S . t . cos i_s

where S is the solar constant (1355 W/m2) and i_s the solar zenith angle. The solar zenith angle is a function of longitude, latitude, time of the day and day number. The daily average global radiation I_g is obtained by integrating the function cos i_sfrom sunrise to sunset. When a pixel is flagged cloudy, the radiation transmission through the cloud (t_c) is estimated from the cloud albedo using a relation from Kubelka-Munk theory. t_c is then used in stead of t in the previous equation.

Net radiation

The net radiation is calculated with:

I_n = (1-A).I_g + L_n

L_n is the net thermal radiation flux, consisting of an upward component emitted by the surface and a downward component emitted by the atmosphere.

L_n= e₀ e_a sT_a4– e₀sT_o4

e₀ is the surface emissivity. A value of 0.9 is assumed. The atmospheric emissivity e_a is estimated with the Brunt equation on the basis of climatic values of air humidity. Below clouds the upward and downward fluxes almost cancel (L_n » 0).

Climatic net radiation

It is convenient to express the long wave radiation flux as:

L_n » e₀ (1- e_a)s T_a⁴ + 4 e s T³ (T₀-T_a) = L_nc + H_r

The first part on the right side is the climatic net long wave radiation (L_nc). The second part is the radiative heat flux (H_r). Our net radiation product is the climatic net radiation. H_ris treated as a part of the sensible heat flux.

Sensible heat flux

The sensible heat may be written as:

H = H_c + H_r = C.v_a (T₀-T_a) + 4 e s T³ (T₀-T_a)

H_c is the convective sensible heat flux. T₀-T_a is the surface-air temperature difference as measured by satellite at noon (sections 5, 6). C is a convective heat exchange coefficient, which depends on the roughness of the earth surface. For bare land (LE=0) a value of C=1 is used. For vegetated surfaces (LE¹0) the value is allowed to grow slowly in line with the vegetation development, until a maximum value of 2.4. The daily average of the sensible heat flux is obtained assuming constant energy partitioning (constant Bowen ratio). The default output of the system is the convective daily sensible heat flux H_c. This product is most convenient for comparison with convective flux measurements, for example based on eddy correlation or scintillometry.

Actual evapotranspiration

The daily actual evapotranspiration is obtained as the residual term in the energy balance equation.

LE = I_n - H

At daily time scale the soil heat flux can be neglected. If pixels are cloudy, the sensible heat flux cannot be determined. However, as discussed in section 8, the net radiation under clouds is estimated. The actual evapotranspiration is then estimated using the energy partitioning determined at the last cloud free day.

Often, the NDVI or another vegetation index is used for drought early warning or as a crop growth indicator. However, the NDVI lags behind the actual evapotranspiration, which represents the real actual crop growth conditions. For example, in the case of a dry spell the spectral behavior of vegetation in the visible and near infrared part of the spectrum (NDVI) will not change immediately, but the actual evapotranspiration will do so, because when experiencing water stress the stomata of leaves respond immediately and will close partly or fully.

Further reading on energy balance mapping, click here.

4) Rainfall mapping

Rainfall mapping is based on the discrimination of cloud top levels in hourly satellite imagery, using temperature thresholds.

Cloud level	IR counts range	Temperature range	Approximate Height range
Cold	< 45	< 226 K	> 10.8 km
High	45 – 59	226 – 240 K	8.5 – 10.8 km
Medium high	60 – 89	240 – 260 K	5.2 – 8.5 km
Medium low	90 – 119	260 – 279 K	2.2 – 5.2 km
Low	119 – air count	279 – air temp. K	< 2.2 km

For each level the Cloud Duration (CD) is counted hourly during a period of 10 days.

Rain gauge data are obtained daily from the WMO-GTS system. A local regression equation is established for each rainfall station;

R = S a_ij CD_ij + d_j

where a_ij is the regression coefficient for the cloud duration at level i at rainfall station j. d_j is the residual at rainfall station j. The regression is based on the observations in station j and the closest 11 surrounding rainfall stations. The regression coefficients and the residual are subsequently interpolated between the rainfall stations for each pixel. Hereafter the rainfall is determined on a pixel by pixel basis.

5) Drought and desertification monitoring

Two drought and desertification indices are generated: the climatic moisture index (CMI) and the soil moisture index (SMI). The CMI was defined by the United Nations Convention to Combat Desertification (UNCCD) in 1994 as:

CMI = L.R / LE_p

LE_p is the potential evapotranspiration in energy units. This quantity is closely related to the net radiation (LE_p » 0.8 * I_n). It can also be obtained from satellite observed air temperatures using the Thornthwaite equation. The result is almost the same.

The CMI indicates a climatic condition. To characterize the actual drought or desertification status of the ground, the relative evapotranspiration or soil moisture index (SMI) is used.

SMI = LE / LEp» 1.25 LE/I_n

The difference between drought and desertification indices is just a matter of time scale. Desertification indices are at least yearly, drought indices are 10 daily or monthly.

6) Crop yield forecasting

Crop yield forecasting is based on a a newly developed crop growth model, which uses the satellite derived radiation and actual evapotranspiration data as input. The model is based on the observation by Monteith, that all crops have about the same growth rate per unit leaf area and that differences in production can be explained by differences in crop geometry and corresponding light interception. This model was extended to include the effects of drought and photosynthetic efficiency.

The gross photosynthesis (P) is given by:

P = c * e * Ig * C * R

where c is a conversion factor, e is the photosynthetic efficiency of light, I_g the global radiation, C the crop coverage and R the crop relative growth. The crop coverage is estimated from the biomass: C = B/B_M, where B_M is the biomass at full cover. The relative growth (R) is derived from the relative evapotranspiration using the formula of Stewart as presented by Doorenbos and Kassam (1979)

R = (1-k) + k * (LE/LE_p)

Here k is a crop specific constant. The light efficiency is calculated with (Rosema et al. 1998)

e = 1 / (1.25 + const * I_g)

Here the constant is also crop specific. The maintenance respiration (M) is expressed as a function of the relative evapotranspiration and is proportional to the biomass (B):

M = c_resp * B * f(LE/LE_p)

The calculation of R and M is done on a daily basis. Every day net biomass increase is added to the total biomass

B_i+1 = B_i + P - M

The simulated crop biomass is used as an indicator of the final economic yield at any stage of the growth cycle. For yield forecasting a relative approach is usually followed. This implies that the biomass development is simulated two times: (1) for the actual evapotranspiration (->B) and once for potential evapotranspiration (->B_p). Subsequently the relative economic yield is assumed to be equal to the relative biomass production.

Y/Y_p = B/B_p

The actual economic yield (Y in kg/ha) can be obtained after estimating the potential economic yield (Y_p) in the area considered. This value may be obtained on the basis of local information such as historic records of reported yield data.

The relative biomass production (B/B_p) stabilizes already halfway the growing season. For this reason satellite derived values of relative biomass production (relative yield) can be used for crop yield forecasting and early warning.

7) Further reading

Doorenbos, J. and A.H. Kassam (1979) "Yield response to water" FAO Irrigation and Drainage Paper 33, FAO, Rome, 193 pp.

EWBMS final report (2001) "European Energy and Water Balance Monitoring System" project ENV4-CT97-0478, European Commission, Research and Development, Programme in the field of Environment and Climate, Theme 3. 148 pages.

For further reading on rainfall and energy balance mapping , order this report on the EWBMS site.

Kondratyev, K.Y. (1969) "Radiation in the atmosphere", New York, London: academic Press.

Monteith, J.L. (1977) "Climate and the efficiency of crop production in Britain", Philosophical Transactions of the Royal Society, London B, 281, 277-294

Rosema A. (1986a) "Results of the Group Agromet Monitoring Project (GAMP)", ESA-Journal, vol. 10, no. 1, 1986, p 17-41.

Rosema A. (1993) "Using METEOSAT for Operational Evapotranspiration and Biomass Monitoring in the Sahel region", Remote Sens. Envir. 45:1-25 (1993)

Rosema A., R. Roebeling, S. Peter, S. Garrido, H.A.R. de Bruin, M.J.M. Saraber, J.N. Roozekrans, F. Garcia, K. Kok, F. Stolle and M.C. Bronsveld (1994) "Assessment and Monitoring of Desertification in the Mediterranean area" (ASMODE), summary final report", project EV5V-CT91-0029, European Commission, Research and Development (XII), Programme in the field of Environment, Topic IV.3.

Rosema, A., R. Roebeling and D. Kashasha (1996) "Using Meteosat for Water Budget Monitoring and Crop Early Warning", in Agrometeorological Applications for Regional Crop monitoring and Production Assessment", Rijks, Terres and Vossen, eds., pp 161-175, EUR 17735 EN

Rosema, A., J.F.H. Snel, H. Zahn, W.F. Buurmeijer and L.W.A. van Hove (1998) "The relation between laser induced chlorophyll fluorescence and photosynthesis", Remote Sens. Environ. 65:143-154 (1998).