Environmental justice and surface temperature: Income, ethnic, gender, and age inequalities

Severe land surface temperature (LST) significantly impacts residents’ thermal discomfort and can be lifethreatening during warm seasons. Therefore, it is essential to identify the inequalities in LST exposure, i.e. the socioeconomic groups who are overor underexposed to LST. There is a knowledge gap in the literature: there is no previous study which differentiates between national-scale inequalities -i.e. inequalities apparent in all location of a country, and the local-scale inequalities -i.e. inequalities existing in some areas of a country. Employing a semi-parametric geographically weighted regression model to study LST in Dutch residential zones in the summer of 2014, the results of this study show that two national-scale inequalities are significant: overexposure of high-income groups and underexposure of the owners of high-value properties. Additionally, eight local-scale inequalities are identified. Among the latter, ethnic inequalities, overexposure of Non-western and Western immigrants, found to be the most severe and frequent at the local scale. Additionally, females are often the second most overexposed to LST at the local-scale. To a lesser degree Rental dwellings and Population age 15− 24 are the second most overexposed of a zone’s population. In the end, the results and policy implications


land surface temperature in warm seasons: a growing challenge
"The 2003 summer European heatwave alone caused up to 70,000 excess deaths over 4 months in central and western Europe," according to the European Climate Adaptation Platform (Climate Adapt, 2021). In summer 2019, "soaring temperatures broke records in Germany, France, Britain and the Netherlands"," Reuters reported (Reuters, 2019), and "several historical records at single locations in France, Switzerland, Austria, Germany, the Czech Republic and Spain," observed the World Weather Attribution (World Weather Attribution, 2019). The severe heatwaves trend is set to continue. A projection of human thermal comfort in 2050 showed that a "substantial increase of heat stress abundance is foreseen in all climate scenarios, for both urban and rural areas " (Molenaar, Heusinkveld, & Steeneveld, 2016). The emergence of European heatwaves coincides with the increase of awareness of the notion of environmental justice. In 2018, the European Environment Agency urged that the inequality of exposure to environmental hazards in the EU had not received the attention it deserves by policymakers (European Environment Agency, 2018). The Netherlands has also shown a high level of vulnerability to heatwaves. Dutch Central Bureau for Statistics (CBS) reported that in the week 30 of 2019 the number of deaths was 15 % higher than an average summer week, nearly 400 more (CBS, 2019).
A spatial phenomenon linking heatwave vulnerability and environmental justice is the land surface temperature (LST). LST has a significant impact on thermal comfort in outdoor and indoor spaces during heatwaves and at summertime (Offerle, Grimmond, Fortuniak, Kłysik, & Oke, 2006) and disproportionately affects the residents with different family types, ages, and dwellings (Oreszczyn, Hong, Ridley, Wilkinson, & Warm Front Study Group, 2006). There is a knowledge gap on the environmental justice of exposure to LST: what are the socioeconomic characteristics of the overexposed groups, and whether or not such socioeconomic profile differs across a country? This study seeks answers to these questions. In the next parts, first, a brief review of previous studies, the knowledge gap, and the aim and the approach of this study are introduced. Subsequently, the methods and the data used in the study are described. Ultimately, the results of the statistical analyses are presented and discussed.

Objective and approach of this study
This research's initial standpoint is that the distinction between environmental injustice's spatial scales is essential for enriching the knowledge and policies over the issue. Such distinction helps identify the causes and set up agendas for mitigating environmental justice by national and local scale studies, visions and policies. This study aims to bridge this knowledge gap by putting forward three research questions: a) Are certain socioeconomic groups of Dutch society over-or underexposed to summertime LST? b) What are the socioeconomic groups who are over-or underexposed in all areas of the Netherlands -hereinafter to be referred to as national-scale inequalities? What are the over-or underexposed groups only in certain areas of the country -hereinafter to be referred to as local-scale inequalities? c) What are the most severe local-scale inequalities across Dutch residential zones?
This study aims to study LST exposure of five categories of socioeconomic groups, who previously found to be unequally exposed to environmental hazards: 1 Non-native ethnic groups, given that their location of living could be affected by their specific position in the housing and labour market (Jesdale, Morello-Frosch, & Cushing, 2013;Mitchell & Chakraborty, 2018); 2 Housing tenure groups, given that heatwave mortality can be higher in the low-quality buildings and access to green is likely to be lower in poor neighbourhoods (Rosenthal, Kinney, & Metzger, 2014;Tan & Samsudin, 2017); 3 Age groups, given that different age groups may dwell in different locations of a city or region and thus be exposed to different levels of heat (Madrigano, Ito, Johnson, Kinney, & Matte, 2015;Vargo, Stone, Habeeb, Liu, & Russell, 2016); 4 Income groups, given that access to green areas could differ across different income groups and income disparities between central and fringe areas can give an income dimension to exposure to urban heat islands (Nesbitt, Meitner, Girling, Sheppard, & Lu, 2019;Park & Guldmann, 2020); 5 Gender groups, given that gender disparities in the labour market and location of living can potentially cause different exposure to summer surface heat (Burkart et al., 2011;Rey et al., 2009).
In the next paragraphs, the method used to analyse the data is presented.

Method
The methodology of this study consists of three steps. At the first step, the interdependencies between the independent variables are controlled for and national-and local-scale inequalities are identified. To do so, an ordinary least square linear regression model (OLS) is developed. It is essential to control that the independent variables' multicollinearity level is small enough that estimated coefficients are meaningful. An indicator for measuring multicollinearity is Variance Inflation factor (VIF). It is commonly accepted that if VIF of all independent variables is smaller than 5, the level of multicollinearity is small enough. The OLS model's result is used for calculation and control of VIF amongst the independent variables. Subsequently, a geographically weighted regression (GWR) is developed. The Geographical Variability Test of the GWR model (developed by Nakaya, Fotheringham, Charlton, & Brunsdon, 2009) is used to identify the national-and local-scale inequalities. The formulation of the OLS model is as follow: where y i is the estimation of average summer LST in the zone i, β 0 is the intercept of the OLS model, β k represents the coefficient of the kth independent variable. x ik and ε i are the value of the kth independent variable and random error term in the zone i. The GWR model is formulated as follow: where (μ i , ν i ) represents the longitude and latitude coordinates of the zone i. β k (μ i , ν i ) and β 0 (μ i , ν i ) are the estimated coefficient and intercept of the kth independent variable specific to the zone i. β k (μ i , ν i ) is calculated by the following formulation: Where β (μ, ϑ) is the unbiased estimate of the coefficient value, and W(μ, ϑ) is an adaptive Gaussian spatial weight matrix: W ij represents the weight assigned to the zone j in the GWR model developed for estimating LST at the zone i. d ij is the geodesic distance between the zone i and the zone j. θ denotes the adaptive bandwidth size -i.e. the maximum number of closest zones included in the GWR model of each of the zones. The optimal value of θ is set in the way that the value of the corrected Akaike Information Criterion (AICc) of the GWR model is minimised. By employing the GWR 4.0 tool (developed by Nakaya et al., 2009), a Geographical Variability Test for each of the independent variables in the GWR is executed. The Geographical variability Test of an independent variable is based on comparing AICc of two consecutive GWR models. In the first model, only the impact of the variable in question is considered geographically invariant, and that of all other variables is set to be geographically variant. In the second model, the impact of all independent variables, among them the variable in question, is considered to be geographically variant. If the AICc of the first model is lower than that of the second model, indicated by positive values of "DIFF of Criterion", the variable in question is a geographically invariant predictor, i.e. representing a national-scale inequality. Conversely, if AICc of the second model is lower than that of the first one, indicated by negative values of "DIFF of Criterion", the variable in question is considered geographically variant i.e. a local-scale inequality.
After identifying national-and local-scale inequalities in the second step of the analysis, a semi-parametric geographically weighted regression (SGWR) is developed. The SGWR model simultaneously estimates the impact of geographically invariant and variant independent variables, i.e., national and local-scale inequalities. The SGWR model is formulated as follows: where β m (μ i , ν i ) is the estimated coefficient of the mth local-scale inequality specific to the zone i, and γ n shows the estimated coefficient of the nth national-scale inequalitya single value estimated for all Netherlands zones. Like the GWR model, the spatial weight matrix is adaptive Gaussian, and the bandwidth size is the number of closest zones that minimises the AICc of the SGWR model. Finally, the three performances of the OLS, GWR and SGWR models are compared, and the model with the best performance is used as the basis for the analysis. The models' performance is compared using three measures: adjusted R2: AICc; the randomness of the spatial distribution of the intercept values (assessed by Moran's Index).
At the third step of the analysis, the estimates of the local-scale inequalities are summarised and compared. First, the coefficients are mapped, and the percentages of the significant (p-value<0.05) positive, significant negative, and non-significant coefficients are presented. Secondly, based on comparing the standardised coefficient values of the local-scale inequalities in each zone, the first and second most severe local-scale inequalities are identified and mapped.

Case study area
This study is conducted on the residential zones of the Netherlands, the so-called wijk in Dutch. The zones are one of the most fine-scale administrative boundaries defined and used by the Dutch central bureau for statistics (CBS, 2014) to report socioeconomic statistics annually. The study is conducted on 2400 out of 2835 zones of the Netherlands. The vast majority of the excluded zones are not residential, e.g. water bodies, ports, natural areas, agricultural lands. A part of excluded zones comprises the areas in which either one of the social and housing data was missing, or enough satellite images could not be retrieved.

Dependent variable
The dependent variable of this study is the average summertime land surface temperature (LST). (Summertime refers to June, July and August.) Using MODIS satellite images (Earthdata, 2020), the LST of five time-periods, each period accounting for eight day and nights, are retrieved. In other words, the data of 40 days and nights are used to estimate the average summertime LST. Using four different sources, the average LST of each time-period is measured at four different overpassing local times: MODIS Terra day (10:30 a.m.), MODIS Terra night (10:30 p.m.), MODIS Aqua day (1:30 p.m.), and MODIS Aqua night (1:30 a.m.). Table 1 summarises the time-periods, overpassing local times, and the sources of the satellite images. Fig. 1 illustrates the average LST at different time-periods and overpassing local times across the Netherlands. The figure shows a peak in LST around the third time-period in late July and gradual cooling towards the end of August.
Ultimately, all LST values of the five time-periods and the four overpassing local times are aggregated at the zones' scale. The size of MODIS imagery and that of residential zones are compatible -i.e. most of the rasters overlap with only one zone. The geometries of the MODIS imagery and the residential zones, however, are different. The former is a 1 km x 1 km raster, whereas the latter consists of irregular administrative boundaries. In this respect, some of the rasters are located at the borders of two or more residential zones, which pose a challenge for the rasters' aggregation at the scale of residential zones. To tackle this challenge, every MODIS raster is divided into 400 smaller rasters (50m × 50m), which are subsequently used to calculate average LST at the scale of residential zones. Fig. 2 illustrates the average summertime LST value aggregated at the level residential zone.

Independent variables
This study uses ten independent variables (Table 2). Provided by the Dutch central bureau for statistics (CBS), two of the independent variables, Western immigrants (%) and Non-Western immigrants (%), describe the ethnic composition of the residential zones. According to CBS, an immigrant is a person with at least one none-Dutch parents. CBS definition includes two generations of immigrants. According to the definition, a second-generation immigrant is a person born in the Netherlands with one or two immigrant parents. A Western immigrant is an immigrant originally from a European country -except Turkey, a North American or an Oceanian country, Japan, or Indonesia. A Non-Western immigrant is an immigrant from Asiaexcept Japan and Indonesia, Latin America, Africa, or Turkey (CBS, 2014).
Three variables represent the housing status of the zones. Rental dwelling (%) is the percentage of the dwellings in a zone which are not occupied by their owners. Property value is the average value of dwellings in a zone. (Note that Property value represents the total value of properties and not the value per square meter.) Building age is the median age of entirely or partially residential buildings, weighted by the buildings' gross floor area. The source for the two first variables is CBS (2014). The third variable is calculated based on the Dutch GIS database of buildings-the so-called 3D BAG database (Esri Netherlands, 2016).
Three variables represent the age group of the residents: Population   Table 1 Time-periods, overpassing local times, and the sources of the satellite images. Each satellite image represents the average LST of 8 days or nights. age 65 or older (%), Population age 15− 24 (%) and Population age 14 or younger (%). The age groups represent the youngest and the oldest residents of the zones, who are likely to be the most vulnerable to heatwaves. One variable shows the gender imbalance in the neighbourhoods: Female minus male (%). The last independent variable is the annual disposable income per head: Income per capita. The variables are taken from or calculated based on the data provided by CBS (2014). Table 2 represents the descriptive statistics of the independent variables.

Identification of national-and local-scale disparities
At the first step of the analysis, an OLS model, with LST as the dependent variable and the ten socioeconomic variables as independent variables, is employed. The OLS model outputs show that the associations between LST and all independent variables are significant at pvalue < 0.05 level. The model produces a high R-square valueexplaining more than 50 % of the LST variation. The Variance Inflation Factor (VIF) of the independent variables shows that all VIF values are well below 5, an often-accepted threshold for VIF, indicating that all   independent variables are reasonably independent and represent a unique impact. Subsequently, a GWR model is employed. The optimal adaptive bandwidth for the model is found to be 110. The Geographical Variability Test is used in order to identify the local-and national-scale inequalities. The result shows that the measurement of the test, DIFF of Criterion, is positive in the case of two independent variablesindicating national-scale inequalities: Based on the GWR model results, on the second step, a SGWR model with national-and local-scale independent variables is developed. Table 3 shows the results of the OLS and GWR models.

Results of the SGWR model and its performance compared to the GWR and OLS models
At the second step of the analysis, a SGWR model is developed. The model distinguishes between the independent variables that exert a similar association with LST across all residential zones of the Netherlandsnational-scale inequalities, and those with spatial variant impactslocal-scale inequalities ( Table 4).
The estimates of the SGWR model significantly outperform those of GWR and OLS models. The SGWR model produces the lowest level of AIC and AICc, indicating the most informative model. The SGWR model produces the highest R-square and adjusted R-square, and therefore offers the best explanation of the associations. The SGWR model's residuals are the closest to a random spatial distribution, indicated by a closer-to-zero Moran's I (Table 5).

Estimates of the national-scale inequalities
The results show that the estimates of both national-scale inequalities are significant at the p-value<0.01 level. The results show that higher levels of Income per capita are associated with overexposure to LST. The higher levels of Property value, however, are associated with underexposure to LST. The absolute magnitude of the standardised coefficient of Property value is found to be almost 60 % larger than that of Income per capita. Fig. 3 compares the standardised coefficients of the national-scale inequalities with the significant standardised coefficients (p-val-ue=<0.05) of the local-scale inequalities. The comparison shows that the association between Property value and underexposure to LST is larger than that of almost all local factor. In other words, there is a noticeable national-scale difference between the residents who can afford more expensive properties, and thus enjoy cooler summer surfaces, and those who cannot do so.
The impact of Income per capita in increasing the exposure to LST is often outsized by the impact of ethnicity. In other words, if there is a significant number of immigrants in a zone, a greater level of inequality is likely to be observed amongst the ethnic groups rather than the income groups. Particularly in the case of the Non-Western immigrants, the observed level of inequality outsizes the inequality based on Income per capita with a large margin. The impact of Western immigrants' presence on increasing the inequality, if significant, is comparable to that of Income per capita. The impact of Income per capita on increasing inequality Table 3 The results of the OLS and GWR models.  significantly outsizes the inequalities amongst the age, genders, and housing tenure groups (Fig. 3). These zones are mainly located at the most urbanised parts of the country, i.e. the so-called Randstad region in the West and some Northern and central zones (Fig. 4a). Non-Western immigrants are found to be overexposed in more than 90 % of the zones (Fig. 4b). in case of the Rental dwellings, in more than 77 % of the zones, no significant disparity is observed. In 20 % of the zones, mainly located in the East of the country, Rental dwellings are associated with LST overexposure (Fig. 4c). Similarly, Building age is not an indicator of inequality in 88 % of the zones. The factor is associated with overexposure only in some of the zones in or at the vicinity of the city of Amsterdam (Fig. 4d).

Estimates of the local-scale inequalities
Population age 65 or older is the only part of the population underexposed to LST in a substantial portion of the zones, 40 %. The zones with significant underexposure consist of the country's Eastern and central zones with high natural areas coverage (Fig. 4e). Population age 15− 24 is overexposed to LST in more than 44 % of the zones. Such zones are mainly located in the less central and northern parts of the Netherlands (Fig. 4f). Population age 14 or younger is not particularly over-or underexposed in most of the zones, 73 %. In the rest of the zones, a mixture of over-and underexposure is observed (Fig. 4g). A significant overexposure of female is observed in 44 % of the zones, indicated by significant positive coefficients of Female minus male. These zones are spread across the Netherlands and comprise urbanised areas such as Amsterdam and the areas in the provinces of Zeeland, Drenthe, Gelderland, Noord Brabant, Overijssel (Fig. 4h). Fig. 5 summarises the percentage of residential zones with over-and underexposed socioeconomic groups at the local scale. Non-Western immigrants are the most overexposed social groups, with LST overexposure in more than 90 % of zones. Three socioeconomic groups are overexposed to LST in 45 %-50 % of the zones: Western immigrants, Population age 15− 24, Female minus male. Three socioeconomic groups are not particularly over-or underexposed in most of the zones: Rental dwellings, Building age, Population age 14 or younger. At last, Population age 65 or older is the only socioeconomic group underexposed to LST in more than 40 % of the zones (Fig. 5).

Estimates of the most severe local-scale inequalities
Comparing the local-scale standardised coefficients, the first and second most severe local-scale inequalities are identified. The most severe local-scale inequality occurs between ethnic groups. In the two most populated cities of the Netherlands, Amsterdam and Rotterdam, and two provinces of Noord-Holland and Friesland, the most severe local-scale inequality is between Western immigrants and the rest of the population. In almost all the other zones, the most severe inequality is between Non-Western immigrants and the rest of the population. In the case of the second, most severe inequality, however, the picture is more diverse. A variety of the factors are identified as the second-largest inequality between the residents, among them gender ethnicity, housing tenure and age group inequalities (Fig. 6). Fig. 7 reports the frequency of the first and second most sever localscale inequalities across the residential zones of the Netherlands. The results show that in 97 % of the zones, the most severe inequalities are observed across the ethnic groups. In almost 70 % of the zones, Non-Western immigrants are the most unequally overexposed part of the population. In 28 % of the zones, the most significant inequality is overexposure of Western immigrants. In the case of the second most severe inequality, Female minus male, as an indicator of the gender imbalance in exposure to LST, is the most common across the zones. To a lesser degree Rental dwellings and Population age 15− 24 are the second most overexposed of a zone's population.  Fig. 3. Estimated standardised coefficient of national-scale inequalities (dotted lines) compared to the significant (p-value < 0,05) local-scale inequalities (box plot).

The national-scale inequalities
This study shows that the inequalities between two types of socioeconomic groups are eminent at the national-scale: Income groups and Property value groups. Presumably, the two inequalities reflect the spatial structure of the Netherlands as a country, the so-called suburbanisation of poverty, and a presence of green as a general characteristic contributing to higher levels of property value in Dutch residential zones. In the following paragraphs, these issues are discussed in detail.
The national-scale inequality between Income groups, i.e. residents with higher incomes being more likely to be overexposed to LST, presumably reflects the gaps between the LST in urban and rural areas and the so-called suburbanisation of poverty. On the one hand, LST is higher in the urban areas, due to the concentration of impervious surfaces, or the so-called urban heat island effect. Such a phenomenon is previously observed in many metropolitan areas among them Paris (Sarkar & De Ridder, 2011), New York (Gedzelman et al., 2003), London (Kolokotroni & Giridharan, 2008). On the other hand, the income gap between core and peripheral areas of many major metropolitan areas is growing, the so-called "suburbanisation of poverty" . Suburbanisation of poverty is also a prominent phenomenon in the two largest Dutch cities of Amsterdam and Rotterdam, where between 2004 and 2013 a significant movement of low-and middle-income occupants from the cities to the urban region is observed (Hochstenbach & Musterd, 2018). In addition to the displacement of the low-income occupants, the income gap between Amsterdam and Rotterdam's central and peripheral areas is partially due to higher levels of social mobility in the central areas (Hochstenbach & Van Gent, 2015). An analysis on the four major Dutch cities of Amsterdam, Rotterdam, Utrecht and The Hague between 1999 and 2014 shows that the growth of average income has been between 5 %-9 % greater in the city cores rather than the metropolitan areas (Modai-Snir & Van Ham, 2020). To sum up, as the cities are gradually gentrifying and poverty is moving to suburbs, high-income inhabitants eventually become more exposed to summer surface temperature than low-income inhabitants.   The national-scale overexposure of high-income residents also reflects the spatial structure of the Netherlands. The west side of the country is "portrayed as one of Europe's most pronounced polycentric mega-city regions," in words of Lambregts, Kloosterman, van der Werff, Röling, and Kapoen (2006, pp. 30). The Randstad region is the most urbanised part of the Netherlands, comprised of the four main Dutch cities of the Amsterdam, Rotterdam, Utrecht and The Hague, and less urbanised areas in the middle, the so-called Green Heart. Presumably, the national-scale inequality between high-and low-income groups, also reflects the contrast between the Randstad and the rest of the Netherlands, in terms of LST and income. On the one hand, the highest LST values during the summer are observed in the Randstad (see Fig. 1), presumably due to the higher presence of impervious surfaces compared to the less urbanised parts of the Netherlands. On the other hand, the average income is significantly higher in the Randstad region's municipalitiesabout 13 % according to the CBS data (own computation based on CBS, 2014). In short, living in a vibrant, mega-city region comes at the cost of exposure to higher LST in summers.
The national-scale inequality characterised by Property value, i.e. residents living in zones with more expensive properties being less exposed to LST, is presumably due to the impact of water bodies, green spaces, and trees on reducing LST and increasing property value. On the one hand, previous studies show that water bodies and green spaces contribute to reducing LST (Klemm et al., 2015;Steeneveld et al., 2014). On the other hand, Luttik, for instance, spotted an "increases in house prices due to environmental factors (up to 28 %) for houses with a garden facing water," particularly if such water bodies are "connected to a sizeable lake" (Luttik, 2000, pp 161). In a study on house price across in the Netherlands, Visser et al. concluded that parks and water bodies' adjacency is associated with higher property values (Visser, Van Dam, & Hooimeijer, 2008). Daams et al. estimated that the value of properties within 250 m from attractive green spaces in Amsterdam could be up to 9 % higher (Daams, Sijtsma, & Veneri, 2019). The effect of attractive natural spaces on property value has been estimated to be between 16 % and 1.6 % in the vicinities of 0.5-7 km (Daams, Sijtsma, & van der Vlist, 2016). A study on property value in the Danish city of Aalborg showed a significant association between green spaces and property value (Panduro & Veie, 2013). Although the impact of open green spaces on property value is vastly affected their size (Cho, Poudyal, & Roberts, 2008), only presence of street trees could appreciate the property value in a neighbourhood (Pandit, Polyakov, Tapsuwan, & Moran, 2013;Siriwardena, Boyle, Holmes, & Wiseman, 2016). In short, the citizens who can afford high-value properties occupy the most adjacent locations to green spaces and water bodies and therefore are less exposed to summer surface temperature.

The local-scale inequalities
A variety of inequalities are in force at the local scale. Non-Western immigrants are the most overexposed group to LST. To a lesser, yet significant extent, Western-immigrants, the population age 15− 24, and females is overexposed to LST. Only One group enjoys underexposure to LST: population age 64 or older. In the following these observations are discussed.
The observed overexposure of non-Western immigrants in about 90 % of the zones is presumably due to particular social and economic processes that determined their living location. There are four main communities of non-Western immigrants: Turkish, Moroccans, Antilleans, and Surinamese. There were two significant inflows of non-Western immigrants to Netherlands: first influx of Turkish and Moroccans as the so-called guest workers in 1960s; second, the influx of Antilleans and Surinamese after independence movements in 1970s (Sleutjes, De Valk, & Ooijevaar, 2018). The non-Western immigrants subsequently dwelled in the country's most urbanised parts, particularly the Randstad region (Musterd & Van Kempen, 2009), typically in neighbourhoods with lowe economic status and affordable housing (Hartog & Zorlu, 2009). Two detailed spatial analyses of the distribution of non-Western population in Amsterdam and Rotterdam's cities showed that the highest concentrations of the population are in central, dense urban areas (Sleutjes et al., 2018;Sleutjes, Ooijevaar, & de Valk, 2019). Non-western immigrants also intend to stay in their current cities and dense neighbourhoods. Boschman and Van Ham, in a study on neighbourhood selection of non-Western immigrants in the Netherlands, found that "ethnic minorities are more likely than others to move to high-ethnic-minority concentration neighbourhoods," and that "own-group effects are also found to be important in explaining neighbourhood selection intend to reside where" (Boschman & Van Ham, 2015, pp. 1155. In short, due to a variety of reasons among them socioeconomic status, own-group effect, and path dependency, non-Western immigrants predominantly live in highly urbanised areas which overexpose them to LST. Western-immigrants are overexposed to LST at almost 50 % of the zones. However, they are less overexposed than non-western immigrants. After the enlargement of the EU in 2004 and 2007, immigration from East to West soared (OECD, 2012). A study by the Social and Cultural Planning Office of the government of the Netherlands showed that most of the immigrants from Poland and Bulgaria, two of the primary origins of western-immigrants, dwelled in the largest cities of the Netherlands (Gijsberts & Lubbers, 2013), and thus are overexposed to LST. Tanis (2020) has observed a similar pattern in European immigrants in German regions: searching for job opportunities, immigrants' first location choice is likely to be a densely populated region. A difference between the Western-immigrants and non-Western immigrants in the Netherlands can explain the lesser overexposure of the former to LST. Western-immigrants appear to be more mobile than non-Western immigrants. In a study titled as "initial and subsequent locations choice of the immigrants to the Netherlands," Zorlu and Mulder observed that immigrants intend to initially dwell in highly-populated, urbanised neighbourhoods, particularly within the four major cities of the Netherlands. Non-Western immigrants are more likely to remain within their initial municipality of dwellings. Western immigrants are almost equally likely to move to another municipality, possibly less urbanised, or stay within the initial one (Zorlu & Mulder, 2008). Presumably, this difference explains the higher level of LST exposure non-Western immigrants compared to western immigrants. Beneath the ethnic layer of inequality, lays a layer of gender inequality. Females are overexposed to LST in almost 45 % of the zones. The unequal gender distribution across urban and rural areas can explain this difference. The socioeconomic data of the Dutch residential zones show that in the zones designated as "very highly urbanised", "highly urbanised" or "moderately urbanised", females outnumber males by respectively 1.39 %, 1.93 %, 0.66 %. Reversely, in the zones designated as "barely urbanised" -referring to rural areasor "no urbanisation" -referring to natural, agricultural or industrial areasmale outnumber females by respectively 0.70 % and 2.70 % (own computation based on CBS, 2014). The latter could be due to the overrepresentation of males in the industry and agriculture economic sectors. According to the Gender Data Portal of The World Bank, in 2019 almost 25 % of employed Dutch male worked in the industry sector, compared to only 6 % of female, and almost 3 % of employed male worked in agriculture, compared to 1.2 % of females. However, females outnumber males in the service sectorwhich are more abundant in urbanised areas: in 2019 92.6 % of employed women worked in the service sector, compared to 72.9 % of males (The World Bank Group, 2020).
The above-mentioned overrepresentation of females in the service sector, which presumably cause females' higher presence in urbanised areas, is in line with the so-called feminisation of labour trend. In her seminal book, Hanna Rosin observed that women dominate 12 out of 15 occupations set to grow in the coming decade, among them those in health services (Rosin, 2012). A study by International Monetary Fund (IMF) concludes that women's participation in the European labour is on the rise and the service sector is the most female intensive sector (Christiansen et al., 2016). In a study on employment in European regions, Guisan and Aguayo (2013) found a significant association between the magnitude of the service sector and women's labour participation in a region. In the case of the highly educated occupations, in particular, the gender inequality between labour forces is to a large extent closed (Evertsson et al., 2009), and the digital labour market is considered to be vastly "feminised" (Duffy & Schwartz, 2018). In short, as a side effect of the feminisation of labour, and particularly services, women outnumber men in the most urbanised areas of the Netherlands and subsequently are significantly overexposed to LST.
Young adults, the population aged 15− 24 years old, are exposed to higher than average LST in 44 % of the zones. This presumably is due to overrepresenting young adults in urban areas. The CBS data on Dutch zones show that the population aged 15− 24 years old is more than 20 % overrepresented in the highly urbanised zones compared to none-urbanised zones, respectively accounting for 14.2 % and 11.6 % of the total population in average (own computation based on CBS, 2014). Overexposure of young adults to summer LST is presumably due to youthification of urban areas, i.e., millennials' willingness to live in highly urbanised and central locations. By identifying young adults clusters in 57 metropolitan areas in the U.S., Moos, Filion, Quick, and Walter-Joseph (2019) concluded that youthification is a widespread phenomenon. A twenty-year longitudinal study of six German and U.S. metropolitan areas observed that central areas have become "younger", whereas suburban areas were ageing (Siedentop, Zakrzewski, & Stroms, 2018). Moos (2016) link the youthification of metropolitan areas with the availability of dwellings too small for households with children. Pfeiffer, Pearthree, and Ehlenz, 2019, pp. 433) explained young adults' intention to live in central areas by their "inside/out, constantly connected lifestyle". Lee, Circella, Mokhtarian, and Guhathakurta (2019) found that the "young, pro-urban" residents of California are the most pro-environmental-policies social group, intending not to own a car and show the highest preference for the use the urban amenities. In short, a combination of lifestyle preferences, the housing market, and anti-car, environmental tendencies has boosted the youthification of highly-urbanised areas, which, in return, increased the exposure of young adults to summer LST.
Compared to other socioeconomic groups, the population aged 65 years or older enjoys LST underexposure in 40 % of the zones. This presumably is related to the overrepresentation of the residents of this age group in less urbanised areas. Population age 65 years or older is about 10 % overrepresented in the none-urbanised zones than the highly urbanised zones, respectively accounting for 18.2 % and 16.6 % of the average population (own computation based on CBS, 2014). A study on the Dutch municipalities' demographic trajectory shows that between 1998 and 2011 the percentage of the population aged 65 years or older grew in the less urbanised municipalities and dropped in the largest municipalities, among them Amsterdam and Rotterdam (Vrieling & Melser, 2013). Overrepresentation of elderlies in rural and less-urbanised areas could be explained by two socio-spatial phenomena: out-migration of rural youth, and in-migration of urban elderly. The tendency of migration to more urbanised areas is found to be large among the rural youth. A study on the rural areas of North-Netherlands showed that between 1995 and 2009, the number of residents aged between 15 and 30 declined by more than 2.5 % (Elshof, Haartsen, van Wissen, & Mulder, 2017). A study on out-migration intention of rural youth in the Netherlands and Belgium identifies attending higher education and seeking new employment opportunities as two principal intentions for out-migration of rural youth (Thissen, Fortuijn, Strijker, & Haartsen, 2010). The second reason for the observed overrepresentation of senior citizens in less urbanised areas could be residential mobility of some of the urban elderlies. A study on residential mobility of Dutch elderlies concluded that "the older people are, the more they tend to move away from larger municipalities" (de Jong & Brouwer, 2012, pp. 37). Notably, there is a higher chance that the people who borne or once lived in rural areas migrate back to their origin in a point of their life (Feijten, Hooimeijer, & Mulder, 2008). In short, LST underexposure of senior citizens could be explained by their overrepresentation in rural areas which presumably related to out-migration of rural youth and in-migration of urban elderlies.

Conclusion
The results of this study urge for studying summer land surface temperature from the lens of environmental justice. High surface temperature in the warm season, this study shows, is, in fact, a spatial manifestation of income, gender, ethnic, and inter-generational inequalities. Over and above being an environmental issue, land surface temperature is a socioeconomic challenge, as its spatial pattern is highly coupled with socioeconomic characteristics of Dutch society, among them suburbanisation of poverty, the housing market and gentrification, half-a-century of ethnic segregation and own-group affiliation of minority groups, the feminisation of labour market, the youthification of metropolitan areas, and ageing of rural areas. Such a new perspective on land surface temperature needs to be multiscalar and multilayer. It needs to be multiscalar because some of the inequalities are eminent at the national scale and thus need to be addressed by national-level policies, and some of the inequalities are specific to particular zones and thus need to be addressed by local-or municipal-scale policies. The approach needs to be multilayer because at the local-scale different types of inequalities are mingled with one another, and location-specific policies need to consider different socioeconomic groups.
This study's results urge for adding a new scope to the Dutch National Climate Accord (Government of the Netherlands, 2019), too. Currently, the policy's focus is the reduction of emissions by energy transition. The link between climate change and the energy consumption is doubtlessly crucial. LST overexposure at the local and national scales will affect energy consumption at local and national scales (Mashhoodi, Stead, & van Timmeren, 2020b) in the case of continued global warming. Nevertheless, the Climate Accord needs to broaden its perspective to add environmental justice into its perspective. The accord needs to account for the projections of environmental and demographic changes, and add protection of most vulnerable socioeconomic groups into its scope.
The results also urge for a new mission for spatial planners: environmental justice in the era of global warming. The Dutch Climate Accord paves the way for a crucial involvement of spatial planners by emphasising the so-called neighbourhood-approach, i.e. neighbourhoods-specific plans involving local stakeholder (Government of the Netherlands, 2019). Spatial planners can play a crucial role by adapting multi-stakeholder, climate-oriented approach to urban development (Gandini, Quesada, Prieto, & Garmendia, 2021), participatory land-use planning (Endo et al., 2017) with consideration of residents' mitigation and adaptation behaviours (Kondo, Mabon, Bi, Chen, & Hayabuchi, 2021), and evaluation of urban plans based on their impact on health (Sahay, 2019). In the coming decades, in other words, spatial planning as a profession faces a challenging task: finding local solutions for a global problem.

Declaration of Competing Interest
The authors report no declarations of interest.