Satellite and Ground Measurements for Studying the Urban Heat Island Effect in Cyprus

An urban heat island (UHI) is a phenomenon whereby an urban area experiences elevated air temperatures due to anthropogenic modification of the environment and is usually more evident at night. During heat waves the local effect of an UHI is superimposed on the regional temperature field and as a result heat stress is enhanced. Both the intensity and the spatial structure of the observed thermal contrast of the UHI depend on various parameters, such as the structure of the urban tissue, the population density and its associated heat release, the land use patterns, the vegetation cover, the surface topography and relief etc. In general terms, the UHI is becoming more intense as city sizes increase. Traditional measurements of the near-surface UHI are based on measurements of the air temperature using urban and rural weather stations or air temperature transects. Thermal satellite sensors, which primarily measure the radiance at the top of the atmosphere in the thermal infrared, retrieve the so called land surface temperature (LST) which is the temperature measured at the Earth’s surface and is regarded as its skin temperature. Given that LST is different from the surface air temperature, a distinction is made in remote sensing studies between surface urban heat island (SUHI) and atmospheric heat island (e.g., Nichol, 1996).

Several studies published in the literature have focused on the use of remotely sensed data for studying the urban heat island effect (Dousset & Gourmelon, 2003; Kato & Yamaguchi, 2005; Lo & Quattrochi, 2003; Streutker, 2002; Tran et al., 2006; Xiao et al., 2007; Yuanbo et al., 2007). Other relevant studies are focusing on the validation of satellite LST retrievals with ground measurements or on the inter-comparison of LST products from different sensors (Mostovoy et al., 2005; Nichol et al., 2009; Retalis et al., 2010). The availability of a multitude of data archives (e.g., from MODIS, ASTER and Landsat TM/ETM+ sensors) with long time-series has recently raised the scientific interest in the relevant field. As a result, several studies have been published on the study of the UHI effect for various cities of the world (Hung et al. 2006; Imhoff et al., 2010; Peng et al., 2012).

This Chapter discusses the urban heat island effect in Cyprus based on both multi-temporal satellite and meteorological data. The necessary information of the study area is provided in Section 2. The description and selection of the heat waves and the analysis of the synoptic conditions favouring the development of heat waves are discussed in Section 3. The development of a Neural Network for the correlation of satellite derived land surface temperature (LST) with ground based air surface temperature is examined in Section 4. The analysis of satellite derived LST for studying the temporal evolution of LST and the deviation of LST (anomaly) from the mean values during a heat wave event are presented in Section 5, while Section 6 refers to the calculation of the mean monthly magnitude of urban heat island (UHI) for the period 2002-2008 and for selected heat wave events.

Study area

The island of Cyprus is located in the eastern part of the Mediterranean Basin. The island is situated between latitudes circles 34° and 36° N, and meridians 32° and 35°E. Cyprus has a typical eastern Mediterranean climate: the combined temperature–rainfall regime is characterized by cool-to-mild wet winters and warm-to-hot dry summers (Michaelides et al., 2009). The climatological annual precipitation of Cyprus is around 500mm. The highest precipitation is recorded in the mountainous areas with 1100mm, while in the coastal areas precipitation is limited to 300-350mm.

From a morphological point of view, the island can be divided into five main morphological regions: (a) The mountainous complex of Troodos located at the center of Cyprus; (b) the mountain range of the Pentadaktylos at the northern part; (c) the central plain of Mesaoria located between of these two mountainous ranges; (d) the hilly areas around the mountainous complex of Troodos; and (e) the coastal plains (see Fig. 1). The coastline of Cyprus is characterized by numerous capes and bays. The narrow coastal plains in the north are covered with olive trees and carob trees, while a short distance from the coast, the northern mountain range (Pentadaktylos) is found, which is a limestone formation and peaks to a height of 1024 meters. At the south and the east of the island there are two salt – lakes.

The Troodos mountain range with a peak at 1951 m covers most of the south-western part and the center of the island. This area is covered almost by forests, mainly pine and other forest trees such as cypresses, oaks and cedars. It is estimated that forests cover about 19% of the total area of the island.

Cyprus is divided into six districts: Kyrenia, Famagusta, Larnaca, Limassol, Paphos and Nicosia (Fig. 2). During the last decade (2000-2010), there has been recorded a dramatic urban expansion (see Fig. 3). As it was found from previous studies (Hadjimitsis et al., 2011), there has been an increase of urban areas of more than 100% compared to late 1980’s and a decrease of 20% of rural areas. These results were derived from an analysis of multi-temporal satellite image classification.

media/image1.jpeg

FIGURE 1.

Main morphological regions of Cyprus

media/image2.jpeg

FIGURE 2.

Districts of Cyprus

media/image3.jpeg

FIGURE 3.

Urban areas of Cyprus shown in red

Heat wave events and synoptic patterns

A strong relationship exists between large scale circulation patterns and regional surface variables such as surface pressure, dynamical rainfall, wind and temperature (Tymvios et al., 2007, 2008, 2010a; Xoplaki et al., 2003). As a consequence, synoptic upper air charts at certain levels comprise a valuable tool for the operational weather forecaster to qualitatively predict occurrences of weather phenomena observed on the ground (e.g., heavy rainfall; see Tymvios et al., 2010a). The height pattern at 500 hPa is often used for this purpose. In order to take advantage of these semi-empirical methods and to simplify the statistical processing, stochastic downscaling methods are often applied to the actual weather patterns in order to generate clusters of synoptic cases with similar characteristics. Weather type classifications are simple, discrete characterizations of the atmospheric conditions and they are commonly used in atmospheric sciences. For a review of various classifications, including their applications, refer to Key & Crane (1986), El-Kadi & Smithoson (1992), Hewitson & Crane (1996) and Cannon & Whitfield (2002).

Heat waves have a distinct impact on society through increased mortality, change in the energy consumption profile and the diversification of social behavior. The severity of the heat events may include the local climatological characteristics, the community design and the individual tolerance to heat. Both the frequency of appearance and the intensity of heat waves are increasing in the Mediterranean area (Founda & Giannakopoulos, 2009).

The eastern Mediterranean climate is characterized by the succession of a single rainy season (November to mid-March) and a single longer dry season (mid-March to October). This generalization is modified by the influence of maritime factors, yielding cooler summers and warmer winters in most of the coastal and low-lying areas. Visibility is generally very good. However, during spring and early summer, the atmosphere is quite hazy, with dust transferred by the prevailing south-easterly to southwesterly winds from the Saharan and Arabian deserts, usually associated with the development of desert depressions (Michaelides et al., 1999). The influence of synoptic types on the urban heat island has been investigated by Mihalakakou et al. (2002) who have also adopted a neural network approach.

The definition for a heat wave recommended by the World Meteorological Organization is “when the daily maximum temperature of more than five consecutive days exceeds the maximum temperature normal by 5°C”. Nevertheless, in most countries, the definition of extreme heat events is based on the potential for hot weather conditions to result in an unacceptable level of adverse health effects, including increased mortality. Also, a threshold in maximum temperature is in practical use in many countries.

These periods of abnormally and uncomfortably hot and (usually) humid weather are very common in the eastern Mediterranean during summer and early autumn. Expert examination of the synoptic patterns on upper air charts may reveal the potential for a heat wave event. In this respect, the research presented here attempts to identify height patterns favorable for heat events by using a neural network classification method, namely, Kohonen’s Self Organizing Maps (see Kohonen, 1990).

DATA

As an indication of a possible heat event, the maximum temperature of Nicosia station in Cyprus was chosen. This station is located within the urban area of the city of Nicosia (35°17’N, 33°35’E, 170m, see Fig. 4) and equipped with traditional instrumentation was operational from 1957 until 2001, when it was upgraded to an automatic station. The database used in this study comprises the maximum and minimum temperature records from this station. The maximum monthly temperature measurements are presented in Fig. 5. Also, for the classification of synoptic patterns, the ERA40 reanalysis for the period of 1958 to 2000 (covering roughly the ERA40 time window) were utilized.

The temperatures database was checked for consistency and homogeneity against measurements from nearby stations while the maximum temperatures were also checked for normal distribution fitting.

METHODOLOGY

The maximum daily temperature at Nicosia station was checked against the climatological monthly average value of the period 1961-1990. If the difference was 5°C or more, then the period was characterized as “possible heat event”. If the subsequent days were also positive against this temperature test for more than three days, then the period was considered as heat event. The heat events were checked against the weather classification patterns in order to identify a connection among particular patterns and heat events. The same procedure was adopted for a difference of 3°C, since events with a 5°C difference are rare even during summer. Special care was taken when checking the last and the first day of the month whereby daily maximum temperature values were subtracted from the average climatological value of the two subsequent months.

media/image4.png

FIGURE 4.

Location of ground stations used

media/image5.jpeg

FIGURE 5.

Box and Whiskers plot of the maximum temperatures in Nicosia (1958 – 2000)

Details of the Self Organizing Maps methodology used for the classification have been presented in Michaelides et al. (2010). The 36-Cluster classification adopted also in the present study has been recently demonstrated by Tymvios et al. (2010b).

RESULTS

The distribution of the heat events in consecutive days for 3°C and 5°C difference is illustrated in Fig. 6. It is clearly evidenced that more than 75% of the events last for 3 to 5 days. Most of the identified heat events occur in the transition periods (i.e., Spring and Autumn). This finding is also supported by the findings in Fig. 5, where the larger variation (the area between 25th and 75th percentile) of the average of the maximum temperatures is given for the same periods. With the exception of the periods 12 to 21 July 1978 (10 days) and 2 to 14 July 2000 (13 days), all incidents lasting more than 10 days for this station occurred in October, November, March, April and May.

Clusters 5 and 34 share most of the heat event occurrences. They are both transition period clusters with similar characteristics, exhibiting a distinctive upper level ridge over the eastern Mediterranean and a deep low to the west of this ridge; Cluster 5 belongs to the cold period and Cluster 34 to the warm period. An example of a Cluster 5 member is illustrated in Fig. 7.

When these clusters appear during early Spring and late Autumn, the heat events last from 8 to 15 days, while when they appear just before or after Summer (May and September) they last around 5 days. Summertime appearances of heat events are equally shared between Clusters 12, 19, 24 and 36, all characterized by warm and dry conditions (Michaelides et al., 2010).

media/image6.png

FIGURE 6.

Distribution of the heat events in consecutive days for 3°C and 5°C difference

media/image7.png

FIGURE 7.

The height pattern at 500hPa (Cluster 5) from 24 November 1962

Although it appears that some Clusters are associated with heat events over Cyprus, the connection between heat events and atmospheric circulation at 500hPa did not give definite results that any of these patterns dominate heat event occurrences (as it was possible to demonstrate in previous studies on rainfall and extreme rainfall events). There might be two reasons for this inadequacy. The first is that the window that was chosen for the classification does not include the synoptic patterns that influence the area sufficiently; the second reason is that, although upper air patterns at 500hPa contribute significantly to the evolution of certain surface features (such as dynamical or extreme rainfall), such an association is not so clear for the temperature field. In the search for associations of the temperature fields with synoptic patterns in the Mediterranean, it is important to consider also the lower parts of the atmosphere.

Future research concerning the connection of the weather classification patterns will be focused into a new, much larger window that will include Northern Africa and the Middle East and a combination of classification of patterns at lower levels of the atmosphere (e.g., 850hPa, 700hPa).

Satellite estimates of temperature versus ground measurements

In this Section, a methodology is presented in which the temperature estimates from the MODIS sensor onboard the Terra Satellite is contrasted with ground measurements. The methodology consists of a neural network approach in which measurements on the ground are used as input to the neural network, whereas, the temperature estimate from the satellite is considered as the network’s output.

The neural network methodology adopted has successfully been implemented before in tackling several climatological problems in Cyprus: the prediction of maximum daily total solar irradiance (Kalogirou et al., 2002), the prediction of the daily average solar radiation (Tymvios et al., 2002, 2005a), the modeling of photosynthetic radiation (Tymvios et al., 2005b) and others.

DATA

For the needs of this research, data from MODIS onboard the Terra satellite have been used. More specifically, the level-3 product MOD11A1 (version 5) for the period 2000-2009 was exploited, at a resolution of about 1km by 1km (0.927km). Using the available Land Surface Temperature (LST) fields derived from MODIS, a time series was established corresponding to the position of ground stations. Wan & Dozier (1996) have developed the Generalized Split Window (GSW) algorithm for the retrieval of LST, using the thermal (infrared) channels of MODIS and under different atmospheric conditions (see also, Wan, 1999, 2008). This algorithm retrieves LST on the basis of emissivities in bands 31 and 32 of MODIS. The accuracy in estimating LST was found to be better than 1K, whereas in most cases it was better than 0.5K (Hulley & Hook, 2009; Coll et al., 2009).

The data base for the surface measurements used in this research consist of the hourly recorded temperature at each of the automatic meteorological stations of the network operated by the Cyprus Meteorological Service (see Fig. 4), in the period 2000-2007. Based on these data, the maximum temperature recorded in the time period 1100 – 1300 UTC (local standard time=UTC+2 hours) was considered as the day’s maximum and was subsequently used in the study.

The training of the neural network implemented requires that there are no missing data in the time series used, because the data are used in groups and are therefore inter-dependent. Therefore, the estimated LST (by the neural network implemented) is based on the data of a whole day and missing values result in the rejection of that day. Following quality control based on the above constraint, the number of automatic stations that were subsequently used was reduced to twelve, as shown in Table 1.

METHODOLOGY

Artificial Neural networks (ANN) are small autonomous computational units (algorithms) with inter-connections which, to a large extent, resemble the functioning of natural computational units, namely, the neurons of the human brain. ANN can be trained and learn through repeated examples so that they can reach conclusions and results without human intervention. Since their invention, ANN covered a wide spectrum of research and disciplines and their application has been phenomenal. A few of the numerous examples of ANN applications are mentioned here: medical systems’ automation for the recognition of malignant tumors, control of military equipment and aircraft, estimation of environmental variables, quality verification in production factories, forecasting of financial indices, weather diagnosis and forecasting etc.

For the implementation of the ANN methodology in the present research, the Multi-Layer Perceptron (MLP) was adopted (see Haykin, 1998). The input to this network is the surface temperature recorded at the ground stations and the output is the temperature estimated by the satellite (LST). Fig. 8 displays the MLP implemented.

media/image8.png

FIGURE 8.

The Multi-Layer Perceptron (MLP) network implemented for the prediction of LST

The data from the twelve ground stations and the respective MODIS estimations of LST were used as follows: 60% were used for Training of the network, 20% for Validation and the remaining 20% as an independent set for Testing.

RESULTS

Table 1 displays the errors in the estimation of LST with the neural network, by using the independent set of data. In this table, the maximum, minimum and average errors along with the standard deviation are shown for each of the ground stations. Overall, the performance of the neural network is considered as very satisfactory. However, there are cases where the error is unacceptable and this requires further investigation.

The relation between input data (ground temperature) and satellite estimated temperature LST (target) is shown in Fig. 9. The results are shown for all the data but also separately for the Training, Validation and Testing data sets. For the Training set, the correlation coefficient is R=0.96991, for the Validation set R=0.89692 and for the Testing set R=0.9145, whereas, for all the data R=0.94747. Based on these findings, the performance of the network in predicting LST is considered as satisfactory.

In this research, an attempt has been made to relate ground measurements of temperature with the temperature as it is estimated from MODIS and develop a neural network methodology that can be used in the estimation of ground temperatures by using the satellite imagery. Although the methodology performs sufficiently, overall, it seems that further refinement is needed in order to improve the approach. The adoption of a single network for all the time series of data seems to limit the application of the methodology. For example, the present single neural network developed for each station does not take into account the large seasonal variations in the parameter concerned. It could be more effective if several neural networks are developed based on seasonally grouped data.

media/image9.png

FIGURE 9.

Neural network performance for the Training, Validation, Test and All data sets

Land surface temperature analysis

The MODIS sensor, onboard Terra and Aqua polar satellites, provides one day and one night image under clear sky conditions. MODIS is particularly suitable for the land surface temperature (LST) product due to its global coverage, radiometric resolution and dynamic ranges for a variety of land cover types and high calibration accuracy in multiple thermal bands.

MODIS LST product is based on the generalized split-window (GSW) algorithm (Wan & Dozier, 1996) using as input the MODIS thermal bands 31 and 32. The parameters in the MODIS GSW depend on the satellite zenith view angles, column water vapor and also on the low atmosphere boundary temperature. The band emissivities rely on the classification-based method (Snyder et al., 1998) according to land cover types in the pixel (Monteiro et al., 2007). Temperatures are extracted in Kelvin; accuracy of 1 Kelvin is yielded for materials with known emissivities (Wan, 1999), while a number of studies have also tested the accuracy of the MODIS LST product with favorable results (Wan, 2002; Wan et al., 2004; Coll et al., 2005; Wan, 2008).

The MODIS Aqua product MYD11A1 (V5) and MODIS Terra product MOD11A1 (V5) – Land Surface Temperature and Emissivity Daily L3 Global 1 km Grid SIN were used. Terra and Aqua overpass times for the study area are considered at approximately 1030 and 1330 UTC for day passes, and at approximately 2230 and 0130 UTC for night passes, respectively.

The use of MODIS LST data for examining the temporal evolution and the retrieved temperature anomaly maps for a heat wave event occurred on 24 June 2007 is presented. Moreover, MODIS LST data are used for calculating the urban heat island (UHI) at four urban areas of Cyprus during the extreme heat wave of August 2010.

MODIS LST TEMPORAL EVOLUTION AND TEMPERATURE ANOMALY MAPS

MODIS LST data were initially used for generating mean monthly climatology LST maps for June in the period of 2003-2008. The mean and maximum Aqua day and night LST values for June are presented in Fig. 10 for the period 2003-2008 for two urban areas (Nicosia, Larnaca) and one rural area (Ag. Marina). The curves show that the mean night LST values for the two urban areas are similar, while for the area of Ag. Marina, the temperature levels are 3-4 °C lower. For all sites, a minimum was observed for year 2005. The situation is different though regarding day LST values. The coastal site of Larnaca exhibited the lowest values among the three areas, while Nicosia and Ag. Marina exhibited similar patterns and temperature levels. The overall trend over time for the three areas showed a positive trend.

The intense heat wave event of 24 June 2007 was next examined in order to study the LST behavior during such events since satellite derived LST is controlled by land cover and topographic effect factors. In Fig.11, temperature anomaly maps, in terms of temperature deviation from the long-term monthly mean values (calculated for the period 2003-2008), are presented for the heat wave event under consideration and for both night and day Aqua passes.

The spatial patterns observed in the temperature anomaly maps are complex. It can be observed that day LST anomaly is more intense (up to 15°C) than the night anomaly. Minimum anomaly is located in the area of the mountain range Troodos (central-eastern part), while the southern part of Cyprus presents higher anomaly values than the northern part. The different values of LST increase are attributed to the difference in the emitted radiance from each land type and/or the urban heat island effect.

media/image10.png

FIGURE 10.

Yearly evolution of the mean Aqua MODIS LST for June (2003-2008) as retrieved from Aqua satellite for three different areas in Cyprus

The amplitude of LST anomaly variation between day and night was examined with the land cover types based on the CORINE 2000 land cover map (Fig. 12). It was found that the mean anomaly amplitude was 2.89-4.05°C for artificial surfaces, 2.87-6.01°C for agricultural areas and 2.81-4.63°C for forest and semi natural areas. However, variations were noticed even in the same category. For example, for artificial surfaces the higher amplitude was noticed for airports and the lower for dump sites. For agricultural areas, the higher amplitude was noticed for pastures and the lower for annual crops associated with permanent crops. For forest and semi natural areas, the higher amplitude was noticed for beaches, dunes and sands and the lower for mixed forest.

A close inspection on the Aqua LST image (Fig. 12) acquired on 24 June 2007 (day pass) depicted that the highest LST values are noticed in areas that are recognized as vulnerable to desertification (Fig. 13). In Cyprus, there are two climatic zones that are considered as sensitive to desertification: the semi-arid, which extends over the larger part of the island and the arid sub-humid, which covers the slopes of the Troodos range and the windward side and higher parts of the Kyrenia range (IACO, 2007).

media/image11.png

FIGURE 11.

Land Surface Temperature anomaly map derived from both day (top) and night (bottom) Aqua MODIS passes for the selected heat wave event

media/image12.jpeg

FIGURE 12.

Simplified CORINE 2000 Land Cover map of Cyprus

media/image13.jpeg

FIGURE 13.

Land Surface Temperature map derived from day Aqua pass for the selected heat wave event

Urban heat island analysis

The variation of the UHI magnitude was examined for the four urban areas of Cyprus based on MODIS Aqua images acquired at night-time (at approximately 0130 local time). The selection of the MODIS Aqua data was based on the criterion that the night-time acquired images allow a more precise LST calculation since there is no incoming solar radiation to change the surface radiation balance, while night-time MODIS LST accuracy has been found to be better than day time (Rigo et al., 2006).

media/image14.jpeg

FIGURE 14.

Areas sensitive to desertification, according to the United Nations Convention to Combat Desertification (IACO, 2007)

The magnitude of the UHI was estimated for each of the four test sites both for the mean monthly period 2002-2008 (Fig. 14) and for selective days of high temperature records of August 2010 (Fig 15). The UHI magnitude was calculated by subtracting the LST value from a rural area (as identified from the position of a pre-selected rural meteorological station) from the respective LST values falling within the urban boundary area of each district on a pixel-by-pixel basis (Tomlinson et al., 2010).

Fig.14 presents the mean monthly maximum UHI intensity for the period 2002-2008 for the four urban areas of Cyprus. As noticed, Nicosia, which is located in the centre of Cyprus, is most vulnerable to UHI during the warm period, when the intensity is recorded above four degrees. On the contrary, the other urban areas (Larnaca, Limassol and Paphos), which are close to the coastline, are lesser affected by UHI during the warm period, with intensities recorded around 1.5 to 3.5°C. These areas also demonstrated high UHI intensities during the cold period.

Next, the spatio-temporal variation of the UHI intensity for each of the urban areas was examined for the period 23 July to 28 August 2010, when high air surface temperatures were recorded (Fig. 14). The temporal variation of the maximum UHI intensity was estimated from the available nocturnal Aqua MODIS images for that period. The results revealed that, for most of the cases, the UHI magnitude curves follow a similar trend. Two major peaks were observed, on 31 July and 25 August 2010.

media/image15_w.jpg

FIGURE 15.

Mean monthly maximum UHI magnitude estimated from MODIS Aqua nocturnal images for the period 2002-2008 for Nicosia, Larnaca, Limassol and Paphos

The spatial variation of the UHI magnitude (Fig. 16) was examined for two dates (31 July and 28 August 2010) and was compared to the mean UHI magnitude as calculated for August for the years 2002-2008. The results derived suggest that, in almost all cases, the spatial patterns of the UHI magnitude observed for each urban area are quite similar to each other with a few variations in the magnitude of intensity due to the severity of the heat wave event. The highest intensities were noticed within the areas of the urban fabric.

The maximum intensities of UHI for each urban area were (a) 31 July 2010: 5.2°C (Nicosia), 3.5°C (Larnaca), 1.9°C (Limassol), and 5.0°C (Paphos) and (b) 25 August 2010: 6.9°C (Nicosia), 3.9°C (Larnaca), 3.1°C (Limassol), and 4.2°C (Paphos). Thus, the deviation form the mean monthly UHI intensities calculated for July and August, correspondingly, were of about 0.6°C and 2.7°C for Nicosia, 0.3°C and 0.8°C for Larnaca, -0.4°C and 1.9°C for Limassol and 3.2°C and 2.8°C for Paphos.

media/image16_w.jpg

FIGURE 16.

Temporal variation of maximum UHI intensity for the four urban areas of Cyprus, as derived from the analysis of Aqua nocturnal data for the period 23 July to 28 August 2010

media/image17.png

FIGURE 17.

UHI estimated from MODIS Aqua nocturnal images for (a) 31 July and (b) 28 August 2010, for the four urban areas of Cyprus, separately

 Concluding remarks

This Chapter commented on the use of Earth observation data along with ground meteorological data for the study of the UHI phenomenon in Cyprus. The synoptic conditions favoring the development of heat wave events were discussed. Neural Network analysis was used for classifying synoptic patterns and relate them with heat events. The majority of the heat events have occurred during the transition periods (Spring and Autumn). However, despite the fact that some clusters can be associated with such phenomena the connection between these events and atmospheric circulation at 500hPa did not give clear results.

Furthermore, an attempt was made in order to correlate ground temperature measurements and MODIS LST data. The results have shown that the methodology can perform sufficiently; however, further refinement is needed in order to improve this approach.

Aqua MODIS retrievals of land surface temperature data were used for studying selective heat wave events. The analysis of LST data depicted the regions that are more prone to such events. The spatial variations of the UHI magnitude was also examined for the major cities of Cyprus, during both mean monthly conditions and for selected events, identifying areas that are most vulnerable.

Related Posts

Comments are closed.

© 2024 Civil Engineering - Theme by WPEnjoy · Powered by WordPress