Split Metropolitan area surface temperature assessment with remote sensing method

Land surface temperature (LST) is a very important parameter for a wide variety of Earth processes. The satellite thermal data enable us to obtain information on land surface temperature. The Landsat Mission is the most important and longest satellite mission with the thermal band capability and a sufficient spatial resolution. Landsat thermal channels have been used in this research to determine the LST within the metropolitan area of Split. The results point to the phenomenon of urban heat islands (UHI), as the simultaneous effect of intensive urbanization in the city of Split metropolitan area, and the observed climate changes.


Introduction
An increasing urbanization causes changes in temperature balance in the densely built up urban areas.Surface air and soil temperature are important factors in urban climatology studies [1][2][3].Thermal remote sensing is used in urban areas for the urban heat island (UTO) observation, land cover determination and classification, and as input data for modelling change of surface atmosphere above urban areas [3][4][5].UHIs are parts of urban areas where the air and soil are warmer compared to the surrounding area (Figure 1).They result from the fact that predominant urban materials such as concrete, asphalt or brick, heat up differently than natural materials such as soil, water or vegetation, and can increase temperature in urban areas so that it is for a few degrees higher than that of the surrounding rural areas [6].In the summer period, maximum air temperature is usually reached at about 4 p.m. [6,7].Highest temperature differences have been observed in summer in cities with over 100,000 inhabitants located in lowland areas or valleys [8].Temperature differences in such areas can reach 12 ºC and more [6,7].Urban heat island effects are bigger, and they locally have more painful and significant consequences, compared to the greenhouse effect.UHI are generally caused by [9]: -Loss of natural vegetation and its replacement with vapour impermeable materials that cause reduced evapotranspiration, lower humidity and higher drought in urban areas; -Dominant urban building materials that increase storage of thermal energy; -Parts of urban geometry, e.g.urban canyons as a predominant structure, especially in downtown areas; -Vehicle exhaust fumes, industry, air conditioning equipment [10].
On the other hand, a rural area does not have the UHI appearance because of high evapotranspiration and large shaded areas created by vegetation and water areas.Factors affecting the UHI [11]: 1. Geographical location: climate, topography, rural environmental; 2. Time: day, season; 3. Synoptic situation (UHI limits): wind, clouds; 4. Urban form: materials, geometry, green areas; 5. Urban functions: energy usage, water usage, pollution; 6. Size of urban community: connection between form and function.
In addition to sun radiation, urban heat islands are mostly influenced by physical properties of dominant building materials [12,13] -High thermal energy emissions -concrete as a dominant urban construction material has high capabilities for thermal energy absorption, while aluminium has low absorption properties -Thermal inertia -property of a material to receive and store thermal energy.Metal and sand have a high thermal inertia (fast heating and cooling), while granite and concrete have a low thermal inertia (slow heating and cooling); -Surface roughness; it is important for UHI creation because high buildings create canyon effects, which increases heating surfaces that absorb heat while also blocking the flow of air masses that cool the space; -Albedo or solar energy reflection coefficient, which is different for each material.
The air and land heating and UHI creation also depend on the land usage, and so the temperature warming differs by type of urban zone (residential zone, industrial zone, tourist zone, suburban zone with vegetation and parks).Industrial zones may be by as much as 7 °C warmer than the surrounding areas [12,14].

Figure 1. Urban heat islands [15]
There are two types or two layers of UHI: surface layer and atmospheric layer [3,15] (Figure 2).Surface heat islands have a higher land surface temperature in urban areas compared to temperature in the surrounding areas.Atmospheric heat islands, characterized by warmer air in urban areas compared to rural areas, can be divided into two separate layers: the canopy layer, which extends in urban areas upwards from the ground surface to the top of trees or buildings, and the border layer extending upwards from the canopy layer to as much as one mile above the ground level.Surface and atmospheric heat islands differ greatly from one another.Surface heat islands are present during the day and night.They are most intense in summer with the daily peak intensity of 3 -8 °C and the nightly peak of 8 -12 °C [16].The method most commonly used for surface land temperature identification is indirect measurement by remote sensing in the thermal region or, farther away, in the infrared region [17].
Based on the analysis of mean annual air temperature sets conducted at 26 weather stations in the Republic of Croatia, Bonacci established that the air temperature suddenly increased at 57.7 % stations, and that this trend was first observed in 1988 [18].
Split Metropolitan area surface temperature assessment with remote sensing method The world's major capital cities have been studying their urban space for decades and so they use new materials and increase vegetation cover whenever possible.Thus they are preparing for climate change by adapting and changing their environment [12,20].On the contrary, similar investigations have not been conducted in the City of Split in the scope of recent urban development interventions, during which the quantity of vegetation was actually reduced (especially In the coastal area where the flow of cool sea air -mistral -has been blocked [4,6,9]), which may contribute to the formation of UHI [11,12].
The aim of this paper is to use remote detection methods to explore micro-climate changes in the land surface temperature in the City of Split metropolitan area, and to determine the existence of UHIs, with their possible temporal and seasonal changes.The research results should sensitise competent physical planning officials about the most vulnerable and critical thermal areas, and provide them with measures for the reduction of surface and air temperature.

Study area
The City of Split is one of the oldest and also the second largest city in the Republic of Croatia.Its geographic coordinates are: 43.5°N and 16.5°E.The city is the economic, administrative, transport (with road, sea, and air links) and tourist centre of the eastern part of the Adriatic coast.The City of Split with its suburban communities (Donje Sitno, Gornje Sitno, Kamen, Slatine, Srinjine, Stobreč and Žrnovnica) occupies an area of 80 km 2 (the city of Split area is 22 km 2 ).Its climate is Mediterranean, characterized by warm and dry summers, with an average air temperature ranging from 21.5° to 25.9 °C.Winters are damp and mild winters with an average air temperature ranging from 7.9° to 10.7 °C.The mean annual air temperature is 16.3 °C.January is the coldest month with an average temperature of 7.9 °C, while July is the warmest with an average temperature of 25.9 °C [21].The lowest temperature measured on 23 January 1963 was -9 ºC, while the highest temperature of 38.6 ºC was measured on 5 July 1950.
An annual average of days with precipitation is 122 days (annual precipitation average is 188 mm) [21].Dominant winds in the winter are bora and sirocco, blowing with the frequency of 35 to 55 % per year [22], while mistral is the dominant wind in the summer period.The Adriatic Sea is a natural water reservoir, with the temperature ranging from 10 to 26 °C, and it is also the most important climate regulator in the wider area [23].The scope of this research covers the metropolitan or wider city area, which includes the urban area, suburbs and suburban communities, and the surrounding towns and local communities, usually formed around the city with a high concentration of population (no less than 100.000inhabitants) [24].The UHI study covers the metropolitan area of the city of Split.The air distance between the most distant points of this metropolitan area (Trogir and Omiš) extends to nearly 70 km along the coastline (Figure 3).According to 2011 census, 178,102 persons reside in the city of Split, while approximately 349,314 inhabitants live in the entire metropolitan area [25].The population density is 2244/ km2 (rising to 7499/km 2 in the city centre) [26].Split economy has been growing at the rate of about 6 %, and so the city is becoming the fastest growing economy in the Republic of Croatia.This sudden growth in economic activity is also reflected in the change of climate in this urban area.

Materials and methods
Scientists use various methods of measurement to identify Urban Heat Islands and to estimate the average temperature of air and ground surface: 1. Direct methods such as: a. Point method (in-situ or Euler's method) for measuring air temperature using fixed digital meteorological stations set up at 1.5-2 meters above the ground surface (recommendation of the World Meteorological Organization -WMO) [3], while the changes in air temperature with altitude are usually measured using sensors mounted on high poles [12].An advantage of such measurements lies in high accuracy of data, which enable creation of timelines (with good time resolution), while principal drawbacks are a small space resolution and high price of data collection Tea Duplančić Leder, Nenad Leder, Željko Hećimović b.Measurement of air temperature using moving sensors (Lagrange method) attached to a vehicle or platform, which moves horizontally along the surface.The changes in air temperature with altitude is measured by sensors fixed on hot air balloons or airplanes [12] 2. Indirect methods; the most common one is the faraway detection by which the land surface temperature is measured (Section 3.1).3. Physical modelling of surface air temperature using small scale models that replace the real surface.These models are used to simulate physical processes occurring in nature and to measure air temperature on the physical model [27].
Surface UHIs identified with indirect methods (faraway detection) are usually represented with thermal images.Atmospheric UHIs are identified with direct measurements on meteorological stations or via moving warmth sensors, and are shown on isothermal maps and temperature diagrams [17].
Many types of models, involving different time and space scales, have been developed in scientific literature in order to identify physical processes that lead to changes in air and soil temperature, and to formation of UHIs.Physical processes occurring in urban areas are a combination of complex and diverse phenomena, which interact with one another.Therefore, the models with simplified physical processes are being built in order to minimise complexity of the model itself, of the process, and of the time needed for data processing.UHI models are being explored at different space scales [27,28]: 1. Space scale of the building, known as the energy model of the building, is being developed as a means to estimate the energy used for heating and cooling; 2. Microscale models analyse interaction of the building with its surroundings, which is a basis for developing microclimate models that are most commonly used in architecture and construction; 3. Physical processes, and hence the appearance of UHIs on the scale of big cities (mesoscale models), are most frequently studied in meteorology and urban climatology.These models are based on the method for solving equations relating to the mechanics of fluids [28].
The remote satellite sensing method, applied on the scale of the wider area of Split, is used in this paper.

Thermal remote sensing
The following platforms are most commonly used for urban and rural surface temperature determination [27] The satellites that are used to observe parts of the Earth surface directly measure the top of the atmosphere or TOA [4] in a very narrow shortwave spectral region (0.2-5μm) and in the longwave spectral region (5-100μm) of electromagnetic radiation, which covers approximately 99.5 % of the total thermal energy emitted from the Earth.The land surface temperature calculation procedure is described in Section 4.
Split Metropolitan area surface temperature assessment with remote sensing method

Meteorological data
For the purposes of this study, air temperature, pressure, and relative humidity data, including soil temperature at 2 cm, were measured for selected meteorological situations at the Split-Marjan meteorological station (φ = 45.5 ºN; λ = 16.4 ºE; h = 125 m) (Tables 1 and 2).The obtained data are the property of the Marine Meteorological Centre in Split and the Meteorological and Hydrological Service in Zagreb.Before further processing, the meteorological data were compared to the open data registered at the meteorological station in Kaštela, so as to avoid rough errors.These data were used to determine atmospheric correction parameters for the selected Landsat 7 and 8 scenes.
According to air measurements, it was established that the summer of 2015 was the warmest summer since the beginning of the meteorological measurements.The continuation of these extreme temperature values is expected [34].

Atmospheric correction
A pre-processing was made before the actual satellite images processing.This pre-processing includes elimination of atmospheric scattering from the scene.As it is widely known, the Landsat TM, ETM+ and TIR provide a thermal infrared image of the Earth's surface in high spatial resolution (each pixel of the scene) based on the TOA radiance, which does not include atmospheric correction.The atmospheric transmission and up-welling and down-welling radiance should be known or modelled for a particular area to enable calculation of atmospheric correction parameters.The method for determining atmospheric parameters based on the target temperature, pressure and humidity is used in the paper.The atmospheric correction for the location of Split, and for the selected date and time, was calculated via the model that uses global atmospheric profiles as modelled by the National Centres for Environmental Prediction (NCEP) [36].This method roughly defines the parameter of atmospheric correction.
Once the global atmospheric profile and the surface temperature are modelled, and using the known pressure and relative humidity information, it is possible to convert the incoming spatial radiance into the outgoing spatial radiance: ( where: τ -the atmospheric transmission ε -the emissivity of the surface L T -the radiance of a blackbody kinetic temperature T L u -the upwelling or atmospheric path of the radiance L d -the downwelling or sky radiance L TOA -the space-reaching or TOA radiance measured by the instrument.
Remark: The radiance is expressed in W/m 2 •ster•μm, while the transmission and emissivity are unitless.
Radiance was converted to temperature using the Planck equation.Details of atmospheric correction calculations are shown in [37,38].The omission of atmospheric correction can result in systematic errors with regard to prediction of surface temperature.It should be noted that the surface temperature can be by 5-10 ºC lower if it is calculated without NCEP atmospheric profile models for calculation of atmospheric correction.Surface meteorological parameters for the meteorological station in Split, and atmospheric correction parameters (atmospheric transmission and upwelling and downwelling radiance) for the selected time and dates are shown in Table 1.Air temperatures measured at the meteorological house of the Split -Marjan meteorological station are shown in Table 2, while the temperature of soil is measured at the some station at the depth of 2 cm.Table 2 shows that, for the same dates, the surface temperature is typically lower than the air temperature during the cold periods of the year, while in summer the surface temperature is significantly higher then the air temperature at mid-day (2 p.m.) reaching the values in excess of 50 °C.High temperatures of air and soil, especially those above 36 °C, can negatively affect human health, because normal cooling of human body is hindered at such temperatures.

Conversion of digital numbers to TOA Radiance
The OLI and TIRS band data must be converted to the TOA spectral radiance using radiance rescaling factors provided in the metadata file [37]: where: L λ -the TOA spectral radiance (Watts/(m 2 *srad*μm)) M L -the band-specific multiplicative rescaling factor from the metadata file A L -the band-specific additive rescaling factor from the metadata file Q cal -the quantized and calibrated standard value of a graphical element (digital number DN).

Conversion to TOA Reflectance
OLI band data can also be converted to TOA planetary reflectance using reflectance rescaling coefficients provided in the product metadata file (MTL file).The following equation is used to convert DN values to TOA reflectance for OLI data: (3) where: ρ λ ' -the TOA planetary reflectance, without correction for solar angle M ρ -the band-specific multiplicative rescaling factor from the metadata file A ρ -the band-specific additive rescaling factor from the metadata file Q cal -the roughly pre-processed and calibrated standard value of a graphical element (digital number DN).

TOA Reflection with correction for solar angle
The TOA reflection with correction for solar angle amounts to: where: ρ λ -the TOA planetary reflectance θ SE -the local sun elevation angle.The scene centre sun elevation angle in degrees is provided in the metadata file θ SZ -local solar zenith angle; θ SZ = 90° -θ SE .
For more accurate reflectance calculations, per pixel solar angles can be used instead of the scene centre solar angles.The per-pixel solar zenith angles for scene centre can be obtained in the Landsat image metadata file.

Conversion to At-Satellite brightness temperature
Temperature radiance conversion must be performed by means of the Planck equation [38].TIRS band data must be converted from spectral radiance to brightness temperature using the thermal constants provided in the metadata file:

Unsupervised classification
The Landsat 8 satellite features two thermal channels and both have a spatial resolution of 120 m.According to the Landsat operating manual, the thermal channel 11 has calibration uncertainties, and so only the thermal channel 10 was used in this study.
After the thermal channel processing, temperatures were classified into 15 classes, approximately by 1 ºC, because that number of classes has proven to be the clearest.Clouds were "masked" from the area, if they were situated within the display area.In addition, the unsupervised classification has proven to be more readable compared to supervised classification, because different dates have different temperature ranges.

Results and discussion
Considerable climate change affect many areas: water resources, economy, human health, energy consumption and social factors [41,42].
The surface temperature at the beginning of August 2000 is shown in Figure 4.For the Split metropolitan area, this temperature ranged at 9:42 a.m.UTC from 15.30 to 59.26 ºC.The areas where highest temperatures ranging from 40 to 50 ºC were recorded are: Karepovac garbage dump site (to the east of Split), and the southern part of Sinjsko polje in the north-eastern and south-eastern parts of the Dugopolje industrial zone.Temperatures ranging from 50 to 59° C can also be detected in a relatively small area, about 500 meters in diameter, near the village of Divojevići in the Lećevica municipality.The temperature in this area is inexplicably much higher than that of the surrounding area, which is most likely due to a minor fire or burning grass in a field.The surface temperature was recorded in mid-July 2015 (figure 5), i.e. in the summer period, just like for the previous scene.The surface temperature for Split metropolitan area, ranging from 22.97 to 53.49 ºC, was determined by interpretation of the thermal channel data.Several areas with the temperature much higher compared to the surrounding zones can also be detected in the scene.Four such characteristic areas, representing specific summer scene details, are shown separately in Figure 7.Here the areas with temperature higher than the surrounding areas, i.e. heat islands, are the quarry areas at the western (Plano) and eastern foothills of Kozjak, and at western parts of Mosor (Perun), and also the Karepovac garbage dump site.Heat islands are also present in areas devastated by fires that are temporarily devoid of vegetation (Plano near Trogir and Lokva Rogoznica to the east of Omiš).The southern part of Sinjsko polje exhibits the same temperature it featured in 2000, i.e. about 45 ºC. Figure 5 shows that the areas around roadways, especially those around the motorway, are much warmer than the surrounding areas, and can therefore also be classified as heat islands.
It can be seen in Figure 5 that the northern part of Split peninsula is much warmer than the southern part of the city, which is due to the cooling effect coming from the sea (sea temperature is much lower than the land temperature, which gives rise to Split Metropolitan area surface temperature assessment with remote sensing method the wind called zmorac that blows from the sea to the shore).This figure shows that the colder zone that follows the coastline is greater in the southern part (250-300 m) compared to the northern part (about 100 m) of the Split peninsula.
It is also evident that the southern half of the peninsula is colder than the northern part.Tea Duplančić Leder, Nenad Leder, Željko Hećimović of weather conditions can clearly be seen on winter scenes, as manifested through surface soil temperatures that differ considerably in front of and behind Kozjak, Mosor and Biokovo mountains.In winter scenes, heat islands mostly occur in denuded rocky areas where temperatures are slightly higher when compared to the surrounding areas.
In addition, temperature differences can clearly be seen on the sea surface, as a consequence of fresh water influx (lower temperature) from the rivers Jadro, Žrnovnica and Cetina, and from submarine springs emerging from the seabed.Those phenomena are not presented because the satellite image processing activity has been focused on land, rather than on sea temperature, which significantly affects the inland temperature.

Conclusion
Determination of soil surface temperature in the metropolitan area of Split by satellite detection has revealed that some parts of the urban area have changed considerably in recent years, as manifested through an increase in soil temperature and occurrence of urban heat islands, while soil temperature in rural areas has been changed to a lesser extent.Over the past several decades the urban area of Split has extended to surrounding areas with vegetation, and so the wider urban area has lost a considerable part of green spaces.A significant part of such green areas has been lost through urbanization, and reduced presence of vegetation in green zones, while in the surrounding rural areas green spaces have been lost as a consequence of forest fires (Plano, Lokva Rogoznica) occurring in summer months.The loss of greenery influences evapotranspiration and thermal characteristics of the wider urban area.It can therefore reasonably be assumed that thermal properties of the city of Split have changed considerably over the past several years as a consequence of rapid and unplanned urbanisation and significant climate changes.
A disturbing occurrence of urban thermal islands has been noted, which is particularly evident in summer months when many people stay in the city as it is a renowned tourist destination.
As the remote detection method has enabled registration of soil surface temperatures in the metropolitan area of Split with a spatial resolution of approximately 100 m, it actually enabled identification of UHI on the scale of 100 m (e.g. a rapidly growing industrial zone in Dugopolje or area surrounding the northern port of Split).Main advantages of this method lie in the fact that the data are free of charge and accessible to everyone, while its drawbacks are the eightday time resolution of data that are registered a little before 10 a.m., which does not coincide with the daily temperature maximum.Satellite images show the area of little less than 200 km, and so it can reasonably be stated that this is the greatest advantage of this method compared to other UHI detection methods.
According to literature data [27,28] UHIs can occur at spatial scales of less than 100 m such as for instance at the scale of one building or several buildings with their surroundings (micro-scale).It can therefore be concluded that the satellite detection of soil temperature used in this paper can not be applied for identification of critical areas (UHIs) on the micro-scale (micro-climate of an urban area).To enable a more accurate determination of critical temperature areas in the future, it would be advisable to combine the satellite detection method with all methods described in Section 3, and especially with direct soil temperature measurements at weather stations, which would be distributed at the spatial scale of less than 100 m.Thus a great quantity of data would be obtained, and these data could be used to check the model based on the method for solving equations relating to mechanics of fluids.
According to results presented in this paper, it can generally be concluded that future research should make use of various methods as this would enable a more detailed study of the occurrence of UHIs, all this in order to change and improve the existing situation in space, to define appropriate environmental protection measures, and to set guidelines for new construction.

Figure 4 . 3 °CFigure 5 .Figure 7 .
Figure 4. Landsat 7 satellite scene for the Split metropolitan area processed on 2 August 2000, surface temperature ranged from 15,3 to 59,3 °C area along the Cetina River.As the motorway was not built in 2000, the image shows only local roads which are warmer than the surrounding area and also the Cetina River course where the temperatures are lower.The last row shows the suburban area and industrial zone of Dugopolje, which has been greatly developed over the past ten years, as can be seen on soil surface temperatures that differ both qualitatively and quantitatively from those registered at the onset of the millennium, i.e. in the period before the global warming.Two winder scenes have also been considered to enable more accurate presentation of surface soil temperatures.The Landsat scene registered on 19 February 2001, with temperatures ranging from -3 to 19 °C, is shown in Figure8.Figure9shows the scene registered on 11 November 2015 with temperatures ranging from 8 to 23 °C.The global change in temperature can immediately be discerned from these two scenes.A distinct distribution

Figure 8 .Figure 9 .
Figure 8. Processed Landsat7 satellite scene from 19 February 2001 showing Split metropolitan area; surface temperature ranged from -3,2 to 19,9 ºC Ground-based sensors that may provide unique perspective of some urban features; they can attain a high temporal resolution and avoid corrections due to atmospheric influences.
[7,17]atellites with extensive spatial coverage, but time intervals may be limited.The data are affected by weather and atmosphere [28] 2. Aircraft with high spatial resolution and a more detailed presentation of urban features.The downsides are high cost, irregular coverage, and non-standardized products 3.The first observation of surface thermal islands (based on satellite sensors) was published in 1972[33].Since then, various satellite sensors have been used for registering temperature of urban areas.Today, many authors mostly use Landsat thermal data to calculate surface temperature.Landsat satellite data were selected in this paper for the land surface temperature determination because they have a good spatial resolution and can be used free of charge (ASTER data are neither free nor available).All other satellite missions have better temporal but worse spatial resolution, which is why their use is less suitable for realizing the basic objective of this study.The Landsat satellite mission has limitations with regard to determination of the land surface temperature.The satellite reaches the zone above the studied area every eight days at 9:42 (UTC).The period of day with maximum surface temperature extends from 14 to 18 hours[7,17].
[35]dard Landsat products, distributed by the United States Geological Survey's Earth Resources Observation and Science Centre (USGS EROS) (http://earthexplorer.usgs.gov/),containroughlyprocessed and calibrated (geocoded) multispectral images recorded by the Operational Land Imager (OLI) and Thermal Infrared (TIR) sensors.Two Landsat 8 OLI/TIR scenes, recorded on 12 July and 17 November 2015, and two Landsat ETM+ scenes, recorded on 02 August 2000 and 19 February 2001, were selected for the purposes of this study.The scenes were selected to represent one characteristic winter and summer meteorological episode (assumption is the largest temperature range) at the beginning of the new millennium (start of Landsat 7 mission) and in the year 2015, which has been declared by the World Meteorological Organization[35]as the warmest year since measurements are made.All pre-processing and processing procedures for Landsat scenes were conducted using the ArcGIS 10.3 program.Landsat satellite scenes were resampled by the cubic convolution method, so as to be oriented to the north and georeferenced in the projection coordinate system of the Universal Transverse Mercator Projection (UTM 33N) on the rotational ellipsoid WGS84 (EPSG: 32633) and stored in GeoTIFF format.Landsat 7 and 8 scenes from the USGS EROS archives are available to all users free of charge and without limitation.The final treatment of the obtained raster maps and data was conducted using the Adobe Photoshop CS6 program.