Identification of thawed and frozen soil state in some Siberia regions by multi-temporal Sentinel 1 radar data in 2017-2018 N V Rodionova 0 , Fryazino, Russia , 141190 2019

The paper deals with the identification of thawed/frozen soils in the topsoil layer for three stations in Siberia: Salekhard, Tiksi and Norilsk by using Sentinel 1B C-band radar data for the period of 2017-2018. Determination of the frozen/thawed soil state is carried out in three ways: 1) by multi-temporal radar data on the basis of a significant in 3-5 dB difference in the backscatter coefficient  0 in the transition of freezing/thawing soil state, 2) by finding the threshold value of  0 at which the temperature in the topsoil layer falls below 00C, 3) by texture features. The first method allows determining the period of time during which the process of freezing/thawing of the soil occurs. The second and third methods allow making local maps of frozen/thawed soils. It is shown that for the studied areas the Spearman correlation coefficient between  0 and air temperature for cross - polarization exceeds the correlation coefficient for co-polarization. The graphs of the AFI (air freezing index) for the period of 2012-2018 are constructed based on the archive data of air temperature for the study areas.

 -

resolution. In contrast to IR and microwave radiometers, synthetic aperture radar (SAR) systems potentially provide information concerning soil freezing dynamics under snow cover at high resolution regardless of cloud cover and time of day and night.

Studies of the thawed/frozen soil state by SAR have been carried out since the early 90-s last century by ERS-1 C-band data [ 5 ], [ 6 ]. Research continued with the missions ENVISAT ASAR [ 7 ], ALOS PALSAR L–band [ 8 ], RADARSAT-2 C-band [ 9 ] and others.

The radar signal from the surface mainly depends on the dielectric constant (DC) of soil which is directly related to the water and ice content, also radar signal depends on surface roughness and vegetation cover.

Active sensors discern frozen soil by detecting variations in liquid water soil content. These variations are related to soil DC values [ 10 ]. The average DC value of dry soil is about 2-3, while that of liquid water is 80. DC soil increases proportionally with volumetric liquid water content. Under critical winter conditions, a greater proportion of soil water freezes which leads to a significant decrease of the average value of soil DC; frozen water has DC of ice (3.2), which is similar to that of dry soil [ 11 ]. The opposite process occurs in spring when the ice within soil melts, which results in an increase in radar backscatter by several dB [ 12 ].

Soil properties (i.e., structure, texture, drainage) affect the soil water content and consequently its freezing dynamics. Mineralized soils freeze more quickly and more deeply than organic soils because of their lesser capacity to store water [ 13 ]. The same relation exists between sandy and clay soils.

The aim of this work is to identify the thawed/frozen soil state in topsoil layer and to construct local maps of frozen/thawed soils in the study areas of some Siberia regions: Salekhard, Tiksi and Norilsk, by using Sentinel 1B radar data. 2. Study area and SAR data 2.1 Salekhard Salekhard is the only city in the country located on the Arctic circle. It is located on the territory of the Yamal-Nenets Autonomous district. The coordinates of the Salekhard meteorological station (WMO 23330) 66031.86’N, 660 40.08’ E. The area of the city of Salekhard is a flat undulating plain on the right Bank of the Poluy river, a major right-Bank feeder to the Ob river (figure 1(a)). It is the area of intermittent permafrost. The depth of seasonal soil thawing/freezing varies from 0.9 m to 3 m and depends on soil lithological composition, moisture and vegetation cover character. The forest covers 8% of the area. The grass cover is poor, the area is mostly covered with mosses and lichens. The territory soil texture is dominated by sandy loam, loam, sand-gravel mixture, but no sand. Mean annual temperature is -6.4°C, the absolute minimum is -54°C, the absolute maximum is + 31°C, the maximum depth of the snow cover is 62 cm.

(a) (b)

Figure 1. (a) Salekhard satellite image, (b) air temperature in degrees Celsius for Salekhard.

Figure 1 shows the satellite image of Salekhard and air temperature plot for the period 26.2.2017 – 9.6.2018 by the rp5.ru data. 2.2. Tiksi Tiksi is an urban-type village, the center of Bulun district of the Republic of Sakha-Yakutia. It is located above the Arctic circle, to the East of the Lena mouth on the shore of Tiksi Bay in the Laptev Sea. Coordinates of Tiksi airport weather station (WMO 21824) 71041.85’ N., 1280 54.18 E (figure 2(b)).

Tiksi Bay is surrounded by low gravelly treeless hills. The sea slightly softens the frosts, the mean January temperature is -37.30C, July and August—+7.60C and +7.70C, respectively. The absolute minimum temperature in Tiksi is -50.50C. The highest possible snow depth is 99 cm. April in Tiksi has the lowest average monthly temperature in Russia, which was in 2006 -27.60C. Figure 3 displays a photograph of the village and air temperature plot for the period 7.2.2017–14.6.2018 by the rp5.ru data. Temperatures were taken at time of site survey of the radar S1.

 2.3. Norilsk-Talnakh

Norilsk is located in the North of Krasnoyarsk region, in the South of the Taimyr Peninsula (WMO ID=23078). Talnakh - district of Norilsk, is 25 kilometers North-East of the center of Norilsk. Talnakh is located on the right Bank of the Norilsk river; (Norilsk-on the left Bank) at the spurs of the Putoran plateau. Talnakh mines are the main raw material base of Norilsk Nickel. At a distance of 2.5 km on Southeast from the Talnakh is R32 CALM ground station with coordinates 69о26’01” N, 88028’03” E. The climate is with negative mean annual temperature -9,800C, long (8 months, from November to May) winter with severe frosts and blizzards, short rainy, cold summer.

 2.4. Sentinel 1B radar data

The study used free open Sentinel 1B (S1) C –band radar data IW (interferometric wide swath) mode, VV and VH polarizations with spatial resolution 10 m. S1 images have been processed by program S1Toolbox and later SNAP [ 14 ]. Data pre-processing includes study area fragment selection and radiometric calibration.

Table 1 describes S1 data used: 1) number of processed images, 2) view angle, 3) radar survey time, 4) number of profile pixels for radar backscatter average.

Place name Salekhard Tiksi Norilsk N=39) for VV polarization (weak correlation) and  S =0.69 (p=10-6) for VH polarization. For Norilsk  S =0.57 (p=6*10-4, N=29) for VV polarization and for VH polarization  S =0.82 (p=0). For Talnakh  S =0.58 (p=4*10-4) for VV polarization and for VH polarization  S =0.74 (p=10-6). For all profiles the correlation in cross - polarization exceeds the correlation coefficient in co- polarization. 3. Methods and results of determining frozen/thawed soil state by radar data

 3.1. The determination of FS state by multi-temporal radar data

The depth penetration e/m wave into the soil is determined by the equation p    ( ' )1/ 2 2   '' , where  the length of the wave,    '  i '' -soil dielectric permeability. For Sentinel 1 the wavelength is 5.4 cm. In the case of frozen soil with  '  5.5 and 0.1   ''  0.5 a penetration depth equals 4  p  20 cm.

Identification of thawed/frozen soil state involves knowledge of changes in soil moisture and DC when the temperature decreases. To understand what changes occur with the soil DC, models of phase transitions (freezing/thawing) of water from the liquid to the frozen state (and vice versa) were created. The best known is the semi-empirical dielectric mixed model Dobson et al. [ 15 ], generalized to negative temperatures Zhang et al. [ 16 ]. The paper [ 17 ] shows the relative DC values for three soil types (silty clay, silt loam, sandy loam) obtained from this model depending on the soil temperature from - 300C to 300C for frequencies from C to Ku-band. It is shown that when the soil temperature drops below 00C, there is a sharp decrease in the values of soil DC for all types of soils and all frequencies. Thus, for the C-band, the jump in the soil DC value is about 20 units for silty clay, 15 for silt loam and about 10 for sandy loam. And as shown by [ 3 ], [ 4 ], the jump in backscatter coefficient is 3-5 dB, what is the basis for identification of thawed/frozen soils by radar data.

Using the time series S1, we can determine the start dates of freezing/thawing of FS by significant changes in backscatter coefficient and then stably low/high  0 values. Figure 8 shows graphs of the difference  0 in the absolute value in dB for the adjacent date radar imagery for Salekhard (a), Tiksi (b) and Norilsk (c) where the maxima differences  0 (shown in the figures by arrows) determine the start date of soil freeze/thaw. thaw freeze thaw freeze (1) thaw

freeze (a) (b) (c)

Based on Dubois et al. [ 18 ] equation that defines the relationship between co-polarization  0 with soil parameters (the real part of the soil DC and the surface roughness), it is possible to determine the jump of the DC during soil thaw/freeze. The differences of the real part of the soil DC are for Salekhard up to 14 units during thawing and up to 5 units during freezing, for Tiksi up to 8 units during thawing and up to 4 units during soil freezing, for Norilsk up to 14 units during thawing and up to 17 units during soil freezing. 3.2. Determining of the FS by surface state factor and  0 threshold value Another way to determine the state of FS only by radar data, without using ground measurements, is given in [ 19 ].

The authors [ 19 ] characterize the state of the soil surface through the surface state factor (SSF): where  S0ummer and  W0 inter - the average  0 values in dB in summer and winter,  0 (t) -the current value of  0 , t–time.

SSF=

 1 2  0 (t)  S0ummer  S0ummer  W0 inter

The authors [ 19 ] argue that the seasonal variation of SSF are highly correlated with seasonal variations in soil temperatures, taking negative values when the soil temperature is below 00C. Then, 5-cm topsoil layer can be classified as frozen if the SSF value measured by the radar becomes negative. Figure 9 shows the SSF graphs for Salekhard and Tiksi.

According to the equation (1) it turns out that the intersection of the value of SSF=0 occurs for the same value  0 (t) , but for different times t , i.e. the value  0 is the same for both the freezing point and the thawing point in the 5-cm upper layer. And this value can be determined for each place by knowing the mean  S0ummer and  W0 inter . Let's call this  0 value by  T0hreshold .

For Salekhard, the value  0

Summer was determined for the period 2.6.2017÷18.9.2017, and for VV polarization the value  0 Summer =-19.2 dB. The value

Summer =-13.6 dB, and for VH polarization  0  0

W inter was determined for the period 30.9.2017÷28.5.2018. For VV polarization  0 W inter =-16 dB, for VH polarization  0

W inter =-22.1 dB. Then, by equation (1) we obtain for  T0hreshold =-14.8 for VV polarization and  T0hreshold =-20.6 for VH polarization.

 For Tiksi, the value  0

Summer was determined for the period 19.6.2017÷11.9.2017, and for VV polarization  S0ummer =-14.9 dB, and for VH polarization  S0ummer =-20.95 dB. The value  W0 inter was determined for the period 23.9.2017÷21.5.2018. For VV polarization  0 W inter =-16.3 dB, for VH polarization  0

W inter =-24.6 dB. Then, by equation (1) we obtain for  T0hreshold =-15.6 dB for VV polarization and  T0hreshold =-22.8 dB for VH polarization.

Using the threshold value T0hreshold , let us construct local maps of thawed/frozen soils, for example, in the area of Salekhard airport for the dates 9.5.2017 and 14.6.2017–thawing, and 18.9.2017 and 12.10.2017–freezing. The resulting images, 2.8*3 km in size, are shown in figure 10. Blue colorvalues  0 < T0hreshold (frozen soil), yellow color–values  0 ≥ T0hreshold (thawed soil). Note the significant decrease in the area of frozen soil on the image for 14/06/2017 during thawing and, consequently, increasing the area of frozen soil on the image for 12.10.2017 during freezing.

 3.3. Determining the status of the FS by textural features

The paper [ 20 ] shows the possibility of distinguishing different types of surface in radar images by the use of Haralick texture features [ 21 ]: ‘contrast’, ‘entropy’ and ‘inverse moment’. ‘Contrast' is a weighted square of brightness levels and serves as a natural measure of the spread of brightness. The ‘inverse moment’ is a measure of local similarity. ‘Entropy’ is a measure of disorder in an image. The area with the highest value of ‘contrast’ distinguishes places with height differences (urban areas and so on), areas with a high value of inverse moment-smooth surfaces. Having built a texture image, it can be used the texture feature 'inverse moment’ as an indicator of frozen/thawed soil, because its value is greater for frozen soils. Figure 11 demonstrates texture images after speckle filtering for test site in Salekhard for dates, as in figure 10 (down), in RGB coding: red-'contrast', green-'entropy', blue'inverse moment’.

9.5.2017

VV 14.6.2017 9.5.2017

Salekhard–thawing

VH 14.6.2017 18/9/2017

VV 12/10/2017 18/9/2017 VH 12/10/2017

Figure 11. Textural images of place under study in Salekhard.

One of the segmentation methods on the basis of paired metrics for extracting homogeneous objects in multi-temporal radar images is given in [ 22 ].

4. Soil freezing/thawing index

Air temperature correlates with soil temperature and soil freezing depth [ 23 ]. Accurate prediction of soil freezing depth is an important aspect of various engineering projects. The real depth of soil freezing is influenced by its mechanical composition, thermal properties, soil moisture, climatic conditions (air temperature, wind speed, precipitation, solar radiation). The freezing depth can be estimated using numerical or analytical methods, but the inputs required in the models are not always available or expensive to obtain. If the input data are available, the results of such models could be quite accurate. In the absence of some input data, it is necessary to estimate how much it affected the accuracy of the result. Of the existing models for predicting the depth of soil freezing, the most wellknown is the Stefan model used today, the Berggren model and the Chisholm and Phang model. The paper [ 24 ] considers the accuracy of prediction of soil freezing depth for these three analytical and semi-empirical models. Since none of the models gives accurate results, the authors [ 24 ] proposed a modified empirical model for different types of soils that requires only daily maximum and minimum air temperatures as input. This model is based on the measured values of the depths of frost penetration in 2010-2011 in the state of Michigan and the calculated values of the cumulative freezing index (CFI): for clay soils with a coefficient of determination of 0.94 (the number of measured data points 29), h  3.3787  (CFI )0.5423 for sandy soils with determination coefficient of 0.91 (the number of measurements 129), n where h is the depth of soil freezing in cm, CFI  Tmaiinr  Tmaiarx  /2 , Tmaiinr and Tmaaixr is daily minimum i1 and maximum air temperature.

The authors [ 24 ] obtained a generalized empirical equation for all soil types, taking into account the measurement of the average thermal conductivity of soil samples in Michigan, noting that the use of their empirical equations for other areas requires calibration.

The maximum depth of soil freezing/thawing for every year is determined by the difference between the maximum and minimum values of CFI. The maximum CFI corresponds to the transition from positive to negative temperatures, the minimum CFI corresponds to the transition from negative to positive temperatures. In the first case, the difference between the extremes corresponds to the assessment of the maximum depth of soil thawing, in the second–the maximum depth of soil freezing. This value is called AFI (air freezing index). Figure 12 shows the AFI values characterizing the maximum depth of soil freeze/thaw for Salekhard, Tiksi and Norilsk for 2012-2018. It is incorrect to talk about the trend in linear approximation here since the coefficient of determination R2 has very small values-from 0.001 to 0.2.

 5. Conclusion

The paper deals with the identification of thawed/frozen topsoil layer state near Salekhard, Tiksi and Norilsk by radar Sentinel 1B C-band data VV/VH polarizations for the period 2017-2018. Determination of frozen/thawed soil state is carried out by three ways: 1) by multi-temporal radar data on the basis of a significant 3-5 dB difference in the values of the backscattering coefficient during the transition to freezing/thawing soil state. The numerical values of the  0 differences are 4.6, 3.4, 5.6 dB for soil thawing for VV polarization and 2.9, 3.3, 4.4 dB for VH polarization, respectively, for Salekhard, Tiksi and Norilsk; 2) by calculating the threshold  t0hreshold at which the temperature in the upper soil layer falls below 00C. The values  t0hreshold allowed creating local maps of frozen and thawed soil in the area of test sites; 3) by the Haralick texture features. It is shown that textural feature 'inverse moment' identifies areas with frozen soil.

It is shown that for the test areas of Salekhard, Tiksi and Norilsk the Spearman correlation coefficient between  0 and the air temperature for cross-polarization exceeds the correlation coefficient for the co-polarization.

The obtained AFI values characterizing the maximum depth of soil freezing/thawing did not reveal a significant trend for the period 2012-2018 according to the air temperature archive data for the studied areas.

 [1] Shiklomanov N I 2005 Physical Geography 4 249 - 263 [2] Duguay C R , Zhang T , Leverington D W and Romanovsky V E 2005 Remote sensing of Northern hydrology , geophysical monograph series 163 91 - 117 [3] Khaldoune J , Van Bochove E , Bernie M and Nolin M C 2008 Proc . IGARSS'2008 (Boston . Mass.) III 382- 385 [4] Khaldoune J , Van Bochove E , Bernie M and Nolin M C 2011 Appl . Environment. Soil Sci. Article ID 193237 1 - 16 DOI: 10.1155/ 2011 /193237 [5] Wegműller U 1990 Remote Sens . Environ . 33 123 - 135 [6] Rignot E , Way J B , McDonald K , Viereck L , Willia C , Adams P , Payne C , Wood W and Shi J 1994 Remote Sens . Environ . 49 145 - 154 [7] Park S-E , Bart A , Sabel D , Wagner W , Naeimi V and Yamaguchi Y 2011 Remote Sens . Environ . 115 3457 - 3467 [8] Mironov V L , Muzalevskij K V , Boike Dj , Shveleva А А , Evtiushkin А V and Filatov А V 2013 Izvestiya vuzov . Fizika 1 1 - 3 [9] Rousseau L P , Magagi R , Leconte R , Berg A and Toth B 2009 Proc. IGARSS'2009 (Cape Town , South Africa) 598 - 601 [10] Hallikainen M T , Ulaby F T , Dobson M C , El-Rayes M A and Wu L K 1985 IEEE Trans . GRS 23 25 - 34 [11] URL : www .iki.rssi.ru/books/2003boyarsky.pdf [12] Ulaby F T , Moore R K and Fung A K 1982 Fundamentals and Radiometry (Dedham, Mass , USA: Artech House) [13] Nyberg L , St¨ahli M , Mellander P E and Bishop K H 2002 Hydrol. Process . 15 909 - 926 [14] URL: https://sentinel.esa.int/web/sentinel/toolboxes/sentinel-1. [15] Dobson M C , Ulaby F T , Hallikaine M T and El-Rayes M A 1985 IEEE Trans . GRS 23 35 - 46 [16] Zhang L , Shi J , Zhan Z and Zhao K 2003 Proc. IGARSS ' 2003 (Toulouse, France) 608 - 611 [17] Jagdhuber Th , Stockamp J , Hajnsek I and Ludwig R 2014 Remote Sens . 6 2008 - 2023 [18] Dubois P C , van Zyl J J and Engman T 1995 IEEE Trans . GRS 33 916 - 926 [19] Mironov V L and Muzalevsky K V 2013 Russian Physics Journal 55 899 - 902 [20] Rodionova N V 2007 Issledovanie Zemli iz kosmosa 34 8 - 14 [21] Haralick R M 1979 Proc . IEEE 67 786 - 804 [22] Plotnikov D E , Kolbudaev P A and Bartalev S A 2018 Computer Optics 42 ( 3 ) 447 - 456 DOI: 10.18287/ 2412 -6179-2018-42-3- 447 -456 [23] Brown W G 1964 Canadian geotechnical J. 1 215 - 226 [24] Rajaei P and Baladi G Y 2015 94st Transportation Research Board Annual Meeting (Washington D.C .) 1 - 15