OF SOIL MOISTURE AND PLANT WATER STRESS MODELS BASED ON SATELLITE THERMAL IMAGERY

12Abstract. The paper analyzes the advantages and disadvantages of the most commonly used groups of models of soil moisture and plant water stress based on satellite thermal imagery. We present a simple proof of linking NDTI and CWSI indicators with plants water stress and quantitative justification for the shape of the points cloud on the chart Ts-NDVI.


INTRODUCTION
The basic problem in the realization of the crop water stress or soil moisture monitoring based on data from the meteorological station is their inadequate spatial density.An increase in the number of weather stations is extremely expensive, so in order to obtain spatial estimates of the weather parameters interpolation methods are used.The interpolation is based on the assumption that the interpolated values changes continuously between the meteorological stations and take values in the range specified by the values recorded.This is often incorrect, in particular in the case of precipitation due to the fact that the size of a single rainstorm cell from a few to a dozen kilometers away (Begum, Otung 2009), are much smaller than the average distance between precipitation posts (precipitation in Poland is measured in about 1200 localizations (IMGW-PIB 2014) which give mean distance between them equal about 16 km).In addition, the moisture content of the soil and plant water stress are influenced also by local conditions outflow associated with the micro relief or heterogeneity of soils which are generalized on soil maps.A potential solution to above described problems have the crop water stress and soil moisture models using satellite thermal imagery.Advantage of thermal imaginary is free of charge availability and prospects for its increasing associated with placing in orbit in 2016 by the European Space Agency ESA first satellite series Sentinel 3 equipped among others in thermal sensors.
For the longest time (since 1982) high-resolution thermal imagery records the satellites of Landsat series.Landsat mission sustainability, its continuity, maintaining a similar spectral ranges and free distribution of archival photos by USGS (United States Geological Survey) has ensured that Landsat are the source of most commonly used images in the scientific work.The disadvantage of images supplied by Landsat is their low temporal resolution -their recovery time over the same place is 16 days.
In turn, the best current source of thermal images recorded with a high frequency (resolution of 1 km) for monitoring soil moisture is sensor MODIS in which are equipped the Terra and Aqua satellites.They are synchronized with each other, so that the satellite Terra orbiting the Earth from North to South passes the Equator in the morning and Aqua satellite orbiting over the same area from South to North, passes the equator in the afternoon.In this way the image of the Earth's temperature for day and night is provided, which can be further used for modeling soil moisture, based on the phenomenon of thermal inertia.In addition, data from the MODIS sensor such as the surface temperature, albedo or NDVI index are distributed free of charge in the form of products with a high degree of processing, in aggregates for periods of 8-day, 16-days and one month.This allows to practically eliminate the most burdensome processes associated with the processing of satellite images, such as removing the cloudy areas and replenishment of their data from the nearest previous photos, for which this area is visible.In practice, for Poland area, in the growing season, it is only able to work with aggregated data from periods of 8 days and occasional replenishing that data with data collected in for longer periods of time.
The usefulness of solutions based on thermal images depend on: the accuracy of soil moisture and crop stress estimation by models, images availability and cost of processing.
Thermal remote sensing (wavelength from 3 μm to 1 mm) All bodies at temperatures above absolute zero (0 [K], -273.16[ 0 C]) emit electromagnetic radiation.The spectral radiance L λ [W m -2 sr -1 m -1 ] of heated  At typical temperature for bodies on the ground of 300 [K] or about 27 [ 0 C] the maximum spectrum value obtained from Planck formula (according to the Wien's displacement law: λmax = 2.898 • 10 -3 / T) lies at wavelengths of 10 [μm] (1 μm = 10 -6 m) which are therefore referred to as thermal TIR (thermal infrared radiation).The possibility of registering radiation (spectral radiance L λ ) in that wavelength by the satellite sensors is used to calculate the surface temperature via Planck's formula.
The spectrum of black body radiation is a good approximation of actual body radiation spectrum.Amendments which are describing derogation of this law are complicated and therefore, in most applications, thermal images are taken into account by assigning to the analyzed surfaces of the constant reduction factor known as emissivity.Emissivity by definition for a black body assumes a value of one, for the white body that does not emit or absorbs radiation a value of zero, and for the remaining bodies referred to as the gray body, intermediate values.
A major problem in the use of the thermal images is the process of wave's absorption by the clouds.Often this prevents analysis of large area at the same time and makes it necessary to connect the thermal images from different days (Landsat return period on the same place is 16 days) varying weather conditions.A connection example picture from two different time periods is shown in fig.2, on the right side where clearly visible connecting line runs from the north to the southwest is.Another problem encountered in the analysis of agricultural area water stress on the basis of thermal images, are mixed crop fields with trees.On thermal images trees always have lower thermal surface temperature and may underestimate the effects of crops water stress.It is associated with higher transpiration of trees, which affects deeper root zone coverage and more turbulent air movements above the canopy caused by their high roughness.This effect is particularly noticeable in thermal images of forests.On the other hand the existence of scattered buildings in rural areas mixed with agricultural fields heat up strongly the area of a single pixel (the spatial resolution of a pixel in thermal band of Landsat7 is 60 m) of thermal image which is sometimes classified as agricultural area and as a result the revaluation of water stress.Both effects -higher temperatures of built-up areas and less of forest are shown in the picture (Fig. 3):

REVIEW OF THERMAL IMAGINARY BASED SOIL MOISTURE AND PLANT WATER STRESS MODELS
The models enable estimation of soil moisture based on thermal satellite images can be divided into two groups: a) Models of heat balance using effect of a stronger heating of the areas where there is a shortage of water in the soil due to the limited under these conditions of actual evapotranspiration which is cooling surface process.
b) Models of thermal inertia based on that the wet areas of heat up and cool down more slowly than of dry areas because water has a high heat capacity which makes that daily temperature amplitude is inversely proportional to the moisture content of the soil.
Key features required input data and the accuracy of prediction models with both of these groups are presented in the following statement: a) Heat balance Models: Evaporation of water is one of the processes in which the surface area reflects the energy supplied by the sun.When evaporation is difficult because the soil is dry, solar radiation causes higher heating of the surface.The models utilize information about the surface temperature obtained from the images in the thermal band to determine the degree of water availability in the soil.The starting point for the construction of models of the heat balance is the heat balance equation of the surface vegetation cover as a consequence of the conservation of energy in unit of time (the energy supplied = energy delivered): Components of the above equation gain a positive sign when the energy is delivered to the plant cover and negative when it is given back.In the equation have been omitted words associated with the heat accumulation by vegetation cover and the accumulation of energy in the process of photosynthesis, as they are in most applications neglecting small.The process of photosynthesis is powered usually by less than a few percent supplied by solar net energy flux Rn (Meyers, Hollinger 2004), but in the case of very dense plant communities such as forests, this proportion may be higher ( When are known: the air temperature Ta, vegetation cover temperature Ts (at the time the images were taken by satellite), wind speed and surface roughness allowing for an estimate of the value of aerodynamic resistance r a , the equation of heat balance and the formula for H allow to calculate actual evapotranspiration ETa: Unfortunately, the calculation of aerodynamic resistance r a is usually a difficult task and introduces the biggest mistake in the estimation of the ETa (Petropoulos 2013) since it depends on a surface heterogeneity (Hasager, Jensen 1999).Difficult is also correct determination of the difference Ts -Ta, and unfortunately the local ground temperature measurements little improve the situation, as is it not known how to extrapolate them to the rest of the area.Models of heat balance can be divided into two groups differ in the way of treating the surface -unilamellar and bilayer (Petropoulos 2013).As the name suggests models of single-layer, surface of the plants and the soil surface are not recognizable and have the same surface temperature Ts.In the bilayer surface models, plant temperature and the temperature of the soil surface are different.This is especially important in case of a rare plant cover when the solar radiation reaches the soil, which in view of the slight soil evaporation outside the range of moisture close to field water capacity (Allen et al. 1998), must be very hot.In such cases recorded by the satellite sensor temperature Ts is much higher than the surface temperature of the plant so not including of this effect in models of single layer leads to an underestimation of evapotranspiration.A disadvantage of two-layer model is much higher complexity algorithms (Petropoulos 2013).
In the group of single-layer models, following models have gained the greatest popularity (Petropoulos 2013):

SEBS -Surface Energy Balance System (Su 2002):
These models use a formula for the ETa, derived from the heat balance equation, substituting for it the surface temperature Ts estimated from the images and the measured air temperature Ta from the meteorological stations.In addition, the original model SEBS is performed by iterative procedure to calculate the parameter of Monin-Obukhov length, characterizing the turbulence properties of the air over the surface of the plant.Unfortunately, SEBS model does not specify how to determine the air temperature Ta outside meteorological stations.

SEBAL -Surface Energy Balance Algorithm for Land (Bastiaanssen et al. a1998, Bastiaanssen et al. b1998):
In models of this type the need to measure the air temperature Ta in the meteorological stations was eliminated through auto-calibration model, assuming a linear relationship between surface temperature Ts and the air temperature Ta: Ts -Ta = α Ts + β, and using the temperature of the two groups of pixels in the thermal images -called cold pixels on vegetated areas with high soil moisture, e.g. on wet grassland for which we know a priori that H = 0: (5) and so-called hot pixel areas of the severed possibility of evapotranspiration, e.g. for lacking plant cover and not waterlogged soils for which we know a priori that ET = 0: (6) By substitution is possible to obtain values of constants α and β which allow using the full heat balance equation to calculate the actual evapotranspiration ETa in rest pixels of image.The disadvantage of the model is a subjec-tive method of selecting pixels for auto-calibration and poor results in strongly carved terrain areas.

METRIC -Mapping Evapotranspiration at high Resolution with Internalized Calibration (Allen et al. 2007):
These models differ mainly from models like SEBAL that, instead of substitution of the a priori H = 0, for cold pixels it uses the following equation: (7) where ETa= Kc•Ks•ET 0 is the actual evapotranspiration calculated as a fraction of reference potential evapotranspiration ET 0 according to well known, used for irrigation FAO56 method (Allen et al. 1998).The choice of the so-called cold pixels is limited to the reference crop (alfalfa) and the choice of the so-called hot pixels is confirmed by the use of the water balance model.The introduction of two layers: the surface of the plants and the soil surface instead of one, makes models of this type, do not underestimate the calculated evapotranspiration, but the disadvantage is the need to measure the temperature of the air.

ALEXI -Atmosphere-Land EXchange Inverse (Anderson et al. 1997):
This is a modification TSEB model which consists on adding the sub model of boundary atmosphere layer so that the properties of this layer are associated with changes in time of the surface temperature.

DisALEXI -Disaggregated ALEXI (Kustas et al. 2003, Norman et al. 2003):
In this model, based on images from two sources of high and low resolution and procedure of matching average value of energy flows counted on the basis of high-resolution images to the values read from the images with low resolution, ground temperature measurements and errors associated with interpolation of them are avoid.
Another group of models, referring to the models of heat balance but usually introduced from empirical observations, are models based on the shape of the cloud of points on the chart where one axis is the surface temperature Ts or any of its function and the second axis some type of vegetation index VI (Fig. 4.).The most commonly used vegetation index is NDVI (Normalized Difference Vegetation Index) (Rouse, JR. et al. 1994), which is designed in such a way that it takes values between -1 and 1.The higher the value, the area is more "green", which means that area has a high content of biomass.Negative index values indicate areas that are devoid of vegetation -uncovered soil, water or buildings.Since the pigment in leaves -chlorophyll -has the property that in photosynthesis strongly absorbs the visible radiation of the red channel R and leaves cell structure strongly reflects radiation from the near infrared NIR , NDVI is based on average values of reflectance in the visible band of red R (λ-0.6 μm) and band NIR (λ-0.8 μm): (8) Models of this type are often used not only to estimate evapotranspiration but also to estimate the soil moisture.and can be calculated on the basis of information on temperature Ts for pixel in given place, and the maximum temperature (Tsmax) and minimum (Tsmin) which have respectively pixels of dry and wet areas on the image.Indicators of this type have the advantage over the others that normalization eliminates errors of absolute shift the level of the measured value -in this case the value of Ts.
Because in a situation of water shortage in the soil, the plants are limiting transpiration by closing the stomata in the leaves, the ratio of actual evapotranspiration ETa to the potential evapotranspiration ET 0 which can be written as: (10) is a measure of the availability of water in the soil for the plants.Analyzing extreme situations: a) in conditions of full water availability: ( b) in conditions of total lack of water: (12) if we denote: α = ρ a Cp/r a and assume that α is uniform over all area of thermal image, we get a system of equations: (13) After substitution ( 13) to (10) we get simple proof that NDTI indicator is equal to ETa/ET0 ratio: (14) which strengthen belief that NDTI can be used to remote measure of the soil water availability for plants.Assumption of uniform α for all area of thermal image is source of practical limitations of use NDTI to situations when we can identify pixels of the same crop with extreme dry and wet conditions.Additionally equation (13) shows that difference Tsmax -Tsmin is maximal for high potential evapotranspiration so index NDTI have small relative errors in high vapour pressure deficit conditions.
We show now that NDTI can be used to explanation of triangular shape of the points cloud on the chart Ts-NDVI.According to the evapotranspiration model FAO56 (Allen et al. 1998) ratio of the ETa/ET0 can be simplify as a product of two evapotranspiration reduction functions: (15) where Ks is reduction caused by soil water deficit in root zone and crop coefficient Kc is reduction associated with degree of crop cover development.If for a given crop the maximum and minimum Kc values are known, corresponding to the surfaces on which it is full crop cover and lack of them, the value of Kc for the crop can be written as (linear mixing model) (16) where fv is a fractional vegetation cover.It can be estimated based on the normalized vegetation index (VI) (we can't use the NDVI index directly because it may take negative values and for the full crop cover reaches a maximum value equal to approximately 0.9): (17) where for the factor n different authors give values from n = 1 (Choudhury et al. 1994) to n = 2 (Carlson, Ripley 1997).In further calculations I'll accept compromise that n = 1.5.Substituting the above formula to the formula for Kc, and assuming that the object of interest are field crops including mainly cereals (average for crops (Allen et al. 1998): Kcmax ~ 1.2; Kcmin ~ 0.4) we get: (18) Because NDTI values are equal to the ratio of the ETa/ET0, it can be write: (19) As you can see the NDVI decreases with an increase in Ts.Also the higher is the availability of water for crop roots (Ks), the lower shifted is NDVI(Ts) curve.This shows that the interpretation of the indicators (NDTI and CWSI), could provide an explanation of the shape of a cloud of pixels for a triangular model and allows for its objective parameter setting in case when we can properly identify pixels of the same crop with extreme dry and wet conditions.Because Ks is known function of soil water amount in root zone and soil retention properties it is possible to use NDTI index for soil moisture assessment purpose.

S-SEBI -Simplified Surface Energy Balance Index:
This model is based on a similar concept as the triangular model (Roerink et al. 2000) with that difference that the axis of the vegetation index was replaced by the albedo (surface 'whiteness').Source of triangular shape is increase of surface temperature when vegetation dries during of water lost and albedo increase, with limitation for surface temperature caused by decrease of absorbed radiation according to decrease of vegetation cover and associated albedo increase.As in the triangular model, calculation of soil water stress based on the cloud of points edges is not fully objective.

b) Models of thermal inertia
The driving force behind changes in soil temperatures is heat transfer.It is widely known that increase in soil surface temperature causes with a certain delay heating her subsequent deeper layers and conversely the decrease in surface temperature results in a reverse process.Daily surface temperature changes are short and cause only soil temperature changes the order of tens of centimeters.Annual temperature changes are long-lasting and cause soil temperature changes in the order of meters.The difference in speed adjustment of the soil temperature to changes in its surface temperature is dependent on soil thermal inertia STI defined by the formula: where C s [J kg -1 K -1 ] specific heat capacity of soil, K [W m -1 K -1 ] is the thermal soil conductivity which determines the flow rate of the heat transfer due to temperature difference and ρ s [kg m -3 ] is density of the soil.Because the air has a low density and low thermal conductivity and water has a high density, capacity and thermal conductivity, the soil thermal inertia is approximately proportional to soil moisture.According to this, thermal inertia models are using the information about changing the surface temperature specified on the basis of two thermal imagery, most often from day and night to estimate the soil thermal inertia and soil moisture.The starting point for the construction of models of thermal inertia is the principle of the conservation of energy for the soil heat flux G : where T(z, t) [K] is soil temperature as a function of time t [s] and the depth z [m] in the soil measured from its surface.The heat flux G [W m -2 ] from soil surface to the soil inside is in accordance to the Fourier-Fick law described by the formula: Substituting the formula (22) for G to the first equation gives the differential equation describing the diffusion changes in temperature of the soil due to the flow of heat.A similar diffusion equation for the flux changes in soil can be obtained by differentiating by z formula (21) and substituting formula (22) for G: Because heat energy flux (G) from soil surface cover to the soil inside is strongly associated with diurnal cycle of temperature changes in soil cover, it can be roughly described as a sinusoidal function of time which fades exponentially with the depth of z: (24) where the Gmax is amplitude of the daily G flux changes, ω is the frequency of these changes and φ is the phase shift depending on the depth z.Relationship with depth z is due to the delay with which the soil on the some depth heats up in relation to its surface.Substituting this solution into the equation on changes the flux G in the soil, gets the formula for determining phase shift: φ=±δz+γ and the formula for the coefficient of specifying how quickly the flux G decreases with increasing depth in the soil: After the calculation of the flux G and its integration we get the formula for T. As had to be expected formulas on the G and T are similar and differ only by phase shift and a constant factor: where A is the value of the albedo.CONCLUSIONS 1.The advantages of models based on the use of thermal images: availability of free images for different temporal and spatial resolution; the existence of multiple models of evapotranspiration and crop stress caused by soil water shortage in root zone based solely on remote data.
2. The disadvantages of models based on the use of the thermal images: low frequency of repetition for images with high spatial resolution (in the case of free Landsat over 2 weeks); possibility of analysis is limited only to cloudless time periods; hardly separation of effects sparse settlement, sparse woody vegetation and shadows of forest -agricultural area border; problems with mixing different crops or land use on the same image pixel especially in conditions of high agricultural holdings fragmentation ; differences between crops albedo, roughness and heat transfer resistance significantly decrease accuracy of models without spatial crop identification and separate crop calibration.
3. Taking into account the need to reduce costly ground-based measurements could be assumed that soil moisture monitoring system should be based solely on remote measurements.Among the models based on heat balance this condition conforms only model DisALEXI.The disadvantage of the Dis-ALEXI model is its great complexity.It makes it difficult at an early stage in the research, in particular for bug tracking and the sensitivity analysis.
4. Among the models based on half-empirical correlation between surface temperature Ts and vegetation indicators should be highlighted standardized NDTI indicator.As shown it is equal to the ratio of actual to potential evapotranspiration so it is a measure of plant water stress (unless there are other stresses), which allows also to monitor soil moisture when known are the retention properties of soil.
5. Other promising method of direct remote soil moisture assessment is Apparent Thermal Inertia ATI index.
] w a ve le n g h t L sp e ctra l ra d ia n ce [at to a temperature T [K] hypothetical body with maximal possibility of emission dependent only on the temperature (called black body), defined as the power of the radiation emitted by unit area of this body in unit solid angle [sr] (1 steradian) per unit length of the emitted waves describes the Planck's law: (1) where: λ [m] -wavelength h = 6.62606957 • 10 -34 [J s] -the Planck constant c = 299792458 [m s -1 ] -the speed of light in vacuum k = 1.3806488 • 10 -23 [J K -1 ] -the Boltzmann constant.

Fig. 2 .
Fig.2.Examples of thermal images (shades of red indicate areas of high temperatures and shades of blue low, urban areas masked) -on the left, in the center there is an embankment of the reservoir tailings "KGHM Iron Bridge", on the right, in the center there is a depression cone of groundwater in the vicinity of the Bełchatów opencast mine (own work based on satellite images Landsat7)

Fig. 3 .
Fig.3.Pulawy and surroundings on 07/07/2013 Landsat7 thermal image.Shades of red indicate warm areas (arable and artificial areas) and shades of blue (waters and forests) -cold (own work) m -2 ] -the net radiation that is energy flux density (power) supplied by the sun to the vegetation cover surface G [W m -2 ] -soil heat flux density that is given back by the vegetation cover to soil H [W m -2 ] -sensible heat flux density that is thermal energy given back by the vegetation cover by convection of heated air from the plant cover λETa [W m -2 ] -latent heat flux density that is energy given back by the vegetation cover on the evaporation of water, ETa [kg m -2 s -1 ] -actual evapotranspiration λ=2.448•106 [J kg -1 ] -latent heat of vaporization.
In the group of two-layer models, the most popular are (Petropoulos 2013): TSEB -Two Source Energy Balance (French et al. 2002, Kustas et al. 2003):

Fig. 4 .
Fig.4.The relationship between the surface temperature Ts and vegetation index NDVI (own work -the points are results of measurements in IUNG-PIB research stations Osiny and Pulki 2013.07.07) you can see the amplitude of the daily changes in temperature soil Tmax is inversely proportional to the thermal inertia STI.In practice relation with soil moisture is analyzed for being approximation of thermal inertia STI index of Apparent Thermal Inertia ATI (Price 1977):