Complexities of urban impacts on long-term seasonal trends in a mid-sized arid city

The variability of surface air temperature is of great importance for both society and the environment, is impacted by global warming and local-scale changes. In the arid eastern part of Washington state of USA, substantial urbanization has transformed the Tri-cities into the state’s third-largest urban cluster. This study utilizes a combination of in situ observations and reanalysis datasets to investigate the influence of land use changes on the region’s 2-meter temperature, revealing local effects that compensate for the background global warming. Within the urban fabric, distinct microclimates have emerged due to varying land use, establishing unique relationships between greenness and temperature alterations caused by land use transitions. Notably, our findings demonstrate that the observed compensating signal in the temperature of farmland locations is primarily attributed to local vegetation increases. Through these observations, this research highlights the urban impacts on local climate, offering valuable insights into the complexities surrounding land use and its consequences on the environment.


Introduction
The global surface air temperature has exhibited a steady increase of approximately 0.96°C from 1850 to 2023 [1], primarily attributed to the increase in greenhouse gases resulting from human activities [2].However, understanding long-term temperature trends at regional and local scales becomes more complex due to multiple factors beyond background climate change.Land use and land cover changes, such as deforestation, forest management, and urbanization, along with alterations in local atmospheric forcing caused by clouds, aerosols, and other factors, can modulate local temperature trends [3][4][5].Notably, urban areas frequently experience higher temperatures compared to non-urban areas, a phenomenon commonly known as the urban heat island effect [6][7][8].However, cities in arid regions may exhibit small or even negative urban temperature signals, referred to as urban cool islands [9].The magnitude of the urban warming or cooling signal depends on various factors, including the diurnal cycle, urban morphology, physical characteristics, irrigation practices, urban extent, waste heat release, and regional climate factors [10][11][12][13].Evaluating the relative importance of the nearsurface air temperature signals in urban regions compared to other compensating factors and feedback is challenging due to the scarcity of long-term urban-scale measurements [14].
Distinguishing the individual impacts of greenhouse gas emissions and land use changes on temperature trends poses a significant challenge, as both factors generally contribute to the overall increase in mean surface air temperature.Previous efforts to estimate the effects of urbanization on temperature have often relied on comparing urban and rural observations, yielding varying results based on the methods used to define urban and rural areas, such as population data or satellite imagery of nighttime light [15][16][17][18][19].A previous seminal study compared temperature data from ∼2000 stations in the US with analysis data based on global model simulations that excluded urbanization effects [20].They found that urbanization and other land use changes account for half of the observed decrease in the diurnal temperature range, with an estimated mean surface warming of 0.27°C per century due to land use changes.
While extensive research has explored the climate impacts of urbanization in larger cities, the effects on smaller and mid-sized cities have received comparatively less attention [21][22][23][24][25].It is essential to note that a significant amount of the global population resides in cities with populations below 300,000 [26,27].These smaller cities often undergo rapid growth and urbanization, resulting in significant changes in the local climate [28,29].Understanding the impacts of urbanization on climate in mid-sized cities is therefore crucial for effective urban planning and management.
This study focuses on the Tri-Cities urban agglomeration located in Washington state, United States.Comprising of three interconnected cities, Richland, Kennewick, and Pasco, the Tri-Cities area is situated at the confluence of the Yakima, Snake, and Columbia Rivers in the semi-arid Columbia Basin of Eastern Washington.Experiencing rapid growth since 2000, the Tri-Cities area has become the state's third-largest city, reaching a population of 300,000 in 2020.Despite this expansion, the Tri-Cities represent a typical understudied city type: a mid-sized city in a semi-arid region with unique characteristics and location.Consequently, this study aims to investigate the influence of urbanization in the semi-arid region of the Tri-Cities on seasonal climate patterns and local warming trends.The subsequent sections describe the datasets used in this study (section 2), analyze the impacts of land use changes on surface air temperature during different seasons (section 3), and provide a summary of the results and implications (section 4).

In-situ observations
Global Historical Climatology Network (GHCN) -Daily is a comprehensive database consisting of climate data compiled from various sources and subjected to rigorous quality review [30].This composite dataset incorporates over 40 meteorological variables, including daily maximum and minimum temperatures, with some station records dating back to the 19th century.The database is regularly updated, and new data is typically available within one to two days of observation.In the Tri-Cities area, Washington state, seven stations are included in the GHCN-Daily dataset, with records spanning from 1950 onwards.Among these stations, three have been in operation for several decades and are located within or in close proximity to the urbanized area, as depicted in figure 1(c).Figure B1 in the appendix B shows the locations of the three stations at a larger spatial scale, while figure B2 depicts the locations of the other five stations.
To visually represent the land use conditions surrounding the three long-operating stations, their locations are shown on Google Maps in figure 1(d).Station 5, situated at the center of the Tri-Cities, is surrounded by a community park and residential area, hence referred to as 'residential.'Station 3, located to the east of the Tri-Cities, is surrounded by commercial and industrial buildings, hereafter referred to as 'concrete.'Station 2 is the furthest from the city among the three selected stations and is characterized by uncultivated arid land and farmland in its vicinity, hence referred to as 'farmland.'These distinct land use settings surrounding the stations will be considered in section 3 to demonstrate the combined response to global warming and local effects stemming from land use changes.

Reanalysis data
The European Centre for Medium-Range Weather Forecast's Reanalysis dataset version 5 (ERA5), developed by the European Centre for Medium-Range Weather Forecasts (ECMWF), is one of the latest global reanalysis datasets available [31].Covering the time period from 1950 to the present, ERA5 assimilates a vast array of ground-based observations, atmospheric sounding data, and remote sensing data to generate a comprehensive set of surface and atmospheric variables.To address land-specific characteristics, ERA5-land was introduced as a land component of the ERA5 climate reanalysis, featuring a higher spatial resolution with a grid spacing of approximately 9 km.This dataset incorporates information on uncertainties for all variables, albeit at reduced spatial and temporal resolutions [32,33].Despite its extensive coverage, ERA5 and ERA5-land lack sufficient in situ measurements from mid-and small-sized cities to represent the effects of land use changes in these areas adequately.Hence, the 2-meter temperature from the ERA5-land dataset is employed as a proxy for capturing large-scale variations in surface air temperature.

Satellite-derived estimates
Vegetation cover plays a crucial role in modulating surface climate, and its significance extends to urban areas [34].In this study, we employ the Normalized Difference Vegetation Index (NDVI), derived from Landsat satellite data, as a reliable proxy for assessing greenness [35].NDVI is a widely used index that effectively captures the presence of live green vegetation at the surface, with documented correlations between NDVI values and the fraction of absorbed photosynthetic active radiation intercepted by plants [36].Furthermore, NDVI has proven valuable in examining the feedback effects of vegetation on local climate dynamics [37].
To estimate greenness, we calculate the mean NDVI for the 30-meter pixel overlaying each station while ensuring cloud-free conditions through rigorous cloud screening of satellite images.Landsat pixels were first removed if any clouds or cloud shadows were detected by the CFMask algorithm before calculating the NDVI [38].This analysis spans multiple seasons and years, specifically from 1985 to 2021.Landsat 5 data is utilized for the period from 1985 to 2012, followed by Landsat 7 for 2013 and 2014, and Landsat 8 for 2015 to 2021, ensuring comprehensive coverage and data continuity [39].

Separating land use changes from other factors
To isolate the specific impacts of urbanization and land use changes from other contributing factors, we adopt a methodology akin to the approach employed by Kalnay et al (2003) [20].This approach involves subtracting the reanalysis datasets from the station datasets, with the annual mean removed from each dataset.This subtraction yields the observation minus reanalysis signal (OMR).
The underlying principle of OMR is that reanalysis data do not incorporate local land surface observations such as for temperature, moisture, and wind speed.As a result, the OMR approach assumes that these variables remain relatively unaffected by urbanization and land use changes and are primarily influenced by assimilated temperature trends [40][41][42][43][44][45][46].It is important to note that the reanalysis data used by Kalnay et al (2003) relied on the NCEP-NCAR dataset, which has been found to exhibit systematic bias in multidecadal variability.This bias can significantly impact the magnitude of OMR signal trends during different time periods [47].
To mitigate this issue, we utilize different reanalysis datasets demonstrating improved agreement with observations by assimilating selected ground-based measurements.It is worth noting that using these alternative reanalysis datasets does not affect the OMR methodology since the ERA5-land reanalysis data employed in this study does not incorporate ground-based measurements from the Tri-Cities area.The urbanization and land use changes in this region occur at a smaller scale than the grid resolution of the ERA5-land datasets, allowing the independence of the OMR analysis from localized ground-based measurements.

Results
The analysis of seasonal averages of the observation minus reanalysis signal (OMR) derived from daily maximum and minimum temperatures reveals a consistent decreasing trend across all seasons and stations with significant correlations (figure 2 and table 1).In contrast, the seasonal averages of the station temperatures themselves do not exhibit significant trends (see figure A1 in the Appendix), suggesting that local effects in the Tri-Cities area are compensating for the overall global warming (figure 1(e)).Notably, the residential station displays a significant and uniform decreasing rate across all seasons, surpassing other stations in terms of the magnitude of this decline.This finding can be attributed to factors such as the expanding area of residential lawns, increased housing construction, and the establishment of community parks in this region, all  contributing to variable surface heating.On the other hand, the concrete stations exhibit the weakest cooling OMR signals, which do not show discernible seasonal variations.This can be attributed to the local warming of concrete surfaces, which offer comparatively less compensation for global warming when compared to areas with vegetation.The farmland station exhibits the highest decreasing rate in summer, with lower rates observed in spring and fall, and no significant trends observed in winter.These seasonal variations in the cooling OMR signals of farmland may be linked to the gradual transition of vegetation type within the farmland region, influenced by the crop growth lifecycle.research has identified several physical mechanisms by which vegetation influences the local climate, such as transpiring water vapor from leaves, increasing atmospheric humidity, and inducing a cooling effect by reducing surface heat loss [34].To assess the contribution of local greenness changes to the observed cooling OMR signals shown in figure 2, we examine the Normalized Difference Vegetation Index (NDVI) derived from Landsat satellite data, which indicates changes in greenness since 1985 (figure 3).In all seasons except for winter, NDVI values are higher in the period 2017-2021 compared to 1985-1989, indicating a general increase in greenness in recent years.The largest differences are observed in summer, while spring and fall show smaller differences.These seasonal differences may be related to the gradual transition of vegetation type in the region from evergreen conifers to more seasonal (deciduous) species, including common lawn grasses.
The averaged NDVI among the nearest ten grids surrounding the three examined stations exhibits variations corresponding to their respective land use types (table 2).In winter, all three locations show a decrease in NDVI, presumably reflecting the prevalence of deciduous lawn grasses that become less visible during this season, contrasting with the region's native coniferous evergreen forests.In spring, summer, and fall, a reduction in NDVI is observed at the concrete locations, possibly due to built-up areas in recent years.Conversely, increases in NDVI are observed at the farmland location during these three warm seasons, which is likely due to the original semi-arid layouts and the crop growth lifecycle.The residential location exhibits the weakest variations in these three warm seasons, likely resulting from a combination of built-up areas (e.g., residential buildings) and grass lawns.
The relationships between NDVI and the OMR signals, as depicted in figure 4, reveal that statistically significant correlations are only observed for the farmland locations.This indicates that in these farmland areas, the variations in vegetation, as represented by NDVI, can explain a significant portion of the observed cooling OMR signals.On the other hand, in concrete locations, where vegetation cover is limited and dominated by built-up areas, the weak cooling OMR signals may arise from other factors that are compensated by the heating effects associated with buildings and infrastructure.On the other hand, the insignificant trends observed in the residential location are not surprising, given the relatively small differences in NDVI values displayed in table 2. This does not imply that land use changes do not contribute to the OMR signals in residential areas but rather suggests that the spatial variations in vegetation and buildings at the 30-meter resolution of the satellite-derived NDVI dataset may not capture the small-scale variations in community parks and residential areas.Therefore,  further investigation is warranted to explore the specific factors contributing to the observed OMR signals in residential areas and understand the complex interactions between land use changes and local climate dynamics.

Discussion and conclusions
Washington State's Tri-cities area has experienced rapid expansion in recent decades, making it the third largest urban agglomeration in the state.Here we study how land use changes have affected the region's 2-meter temperature, and our findings indicate that the area has consistent cooling signals from land use changes.These cooling signals compensate for the changes due to global warming and can partly be attributed to the increased greenness resulting from land use changes.We also show that different microclimates have emerged due to variations in land use changes, resulting in unique relationships between greenness and the temperature changes caused by land use modifications.
The analysis of observation minus reanalysis signal (OMR) derived from daily maximum and minimum temperatures reveals a consistent decreasing trend across all seasons and stations in the Tri-Cities area.The residential station exhibits the most significant and uniform decline, attributed to factors such as expanding residential lawns, increased housing construction, and establishing community parks.Concrete stations display weaker cooling OMR signals, likely due to higher surface sensible fluxes from concrete compared to areas with vegetation.The farmland station shows the highest decreasing rate in summer, with variations across seasons linked to the transition of vegetation types within the farmland region.Normalized Difference Vegetation Index (NDVI) analysis indicates an overall increase in greenness since 1985, particularly in recent years, with larger differences observed in summer.NDVI variations among the nearest 10 grids surrounding the specific stations align with their respective land use types, including a decrease in winter and reductions in spring, summer, and fall in concrete areas, while farmland areas show increases during these warm seasons.The relationships between NDVI and OMR signals demonstrate statistically significant correlations primarily for the farmland locations, indicating that vegetation greenness may play a pivotal role in generating the observed cooling OMR signals.Further modeling-based analysis in the future may be required to confirm and refine this hypothesis and explore any additional factors that may contribute to the negative OMR signals.
The findings highlight the significant influence of vegetation and land use changes on the observed OMR signals.While the farmland areas show clear associations between vegetation changes and cooling OMR signals, the limited vegetation cover and dominant built-up areas in concrete locations contribute to weaker cooling signals.The insignificant correlation between NDVI and OMR signals observed in the residential areas may be attributed to the small-scale spatial variations not adequately captured by the 30-meter resolution of the NDVI dataset.Further investigations are warranted to delve into the specific factors driving the observed OMR signals in residential areas and to gain a deeper understanding of the complex interactions between land use changes and local climate dynamics.
There has been substantial discussion on the role of urbanization on measured temperature trends in many parts of the world [46,48,49].This discussion has primarily focused on whether urbanization, through urban expansion and densification, contributes to long-term warming trends by impacting the footprint of nearby weather stations.Here, we find an interesting counterexample that shows the opposite urban impact.For this semi-arid region, we find that if urbanization were taken into account, it would actually decrease the overall warming trends due to local urban greening over time.These findings are particularly intriguing because this result is consistent with broader definitions of urbanization that are not restricted to mere increases in built-up areas.Urbanization can encompass a wide range of human-mediated land cover modifications, including practices like lawn irrigation and tree planting efforts, each characterized by specific traits unique to their respective cities and regions [9,13,50,51].These modifications may also yield similar observations of cooling due to increased greenery in cities [52][53][54][55][56][57].In general, small or negative urban heat islands (referred to as urban cool islands) have been often observed for cities in arid and semi-arid areas, including a general cooling in arid cities over time [58][59][60].Our findings stress the importance of a nuanced understanding of the relationship between urbanization and rising temperatures and underscore the importance of considering multiple factors when studying climate change.
Also, it is important to acknowledge certain limitations within the scope of this study.Firstly, the method employed in this research relies on statistical techniques to infer the cooling effects of anthropogenic activities, potentially linked to local vegetation increases.Future research could provide more comprehensive details, including sensitivity studies using modeling approaches, which would enhance our understanding of the underlying mechanisms involved and driving factors of this cooling effect [57].Secondly, expanding the focus to include similar cities would bolster our findings' generalizability and contribute to a more holistic understanding of the dynamics at play.Thirdly, urbanization and land use changes may interact non-linearly with global warming [51].Future studies could extend to investigate the urban cooling effect in a warmer climate.These limitations underscore the need for continued research in this domain, with the potential for more robust and nuanced insights in subsequent studies.

Figure 1 .
Figure 1.(a) Population in the Tri-cities area.Two of the three cities are the major area in Benton County, and Pasco is the major area in Franklin County regarding population.(b) Operational period of the stations in the Tri-cities area.(c) Location of the stations.Shading in figure c is the leaf area index from ERA5-land under the constant assumption.Grey lines in (c) donate the rivers.(d) Satellite images from Google Maps with the zoom-in figures of the three long-time operated stations.(e) Time series of annually-mean 2 m temperature from ERA5-Land and PRISM.

Figure 2 .
Figure 2. Long-term trend in the OMR signals for daily maximum air temperature for (a) winter (DJF: December, January, and February), (b) spring (MAM: March, April, and May), (c) summer (JJA: June, July, and August), and (d) fall (SON: September, October, and November).T st is the 2 m air temperature anomalies from stations respective to their 70-years annual cycle.T era5l is the 2 m air temperature from ERA5-land reanalysis datasets after removing their annual cycle.We only show the linear regression lines with p-values smaller than 0.05.The correlation coefficients are shown in table 1.

Figure 3 .
Figure 3.The density of the differences between the averaged NDVI during 2017-2021 and 1985-1989.

Figure 4 .
Figure 4. Correlations between the NDVI and OMR signals.The solid line is the linear regression of the correlations for those with p-values smaller than 0.05.The numbers are the coefficients of the linear regression.

Table 1 .
Correlation coefficients of OMR signals in figure 2. The 'NaN' in the table indicates indicates that there are no significant trends.

Table 2 .
Difference of NDVI between 2017-2021 and 1985-1989 for the three stations.