How the Interaction of Heatwaves and Urban Heat Islands Amplify Urban Warming

Gemechu Fanta Garuma ^1,2,*

1. University of Quebec at Montreal (UQAM), Faculty of sciences, Institute of Environmental Sciences, Canada

2. Entoto Observatory and Research Center, Atmospheric and Climate Science Unit, Department of Space Science and Application Research at the Ethiopian Space Science and Technology Institute, Addis Ababa, Ethiopia

* Correspondence: Gemechu Fanta Garuma

Academic Editor: Zed Rengel

Special Issue: Climate Change and Land

Received: March 29, 2022 | Accepted: May 24, 2022 | Published: June 06, 2022

Adv Environ Eng Res 2022, Volume 3, Issue 2, doi:10.21926/aeer.2202022

Recommended citation: Garuma GF. How the Interaction of Heatwaves and Urban Heat Islands Amplify Urban Warming. Adv Environ Eng Res 2022; 3(2): 022; doi:10.21926/aeer.2202022.

© 2022 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

Urban heat island (UHI) remains one of the top consequences of urbanization, population growth and climate change. It has been studied extensively since it was identified by Luke Howard in the mid 1980’s [1]. Since then both qualitative and quantitative approaches have been followed to study UHI [2] encompassing theoretical, experimental, and numerical techniques. Key synoptic, mesoscale and local scale meteorological factors and the non-linear interactions controlling it have been the subject of many urban climate researches. Statistical models that relate urban heat islands with population growth and fraction of impervious surfaces had mostly been used before the 1990s. However, in the early 1990’s developments of numerical climate models implementing urban surfaces over a wide range of scales and details appeared in literature [3]. The level of detail and complexity required for a particular purpose has been challenging for the urban climate modeling communities. According to Baklanov’s [4] “fitness-for-purpose” approach, for meteorological purposes such as the quantitative studies of the urban heat island (UHI) effects, heatwaves, precipitation and surface runoff, the single layer urban canopy models, SLUCMs are sufficient. But, for more detailed indoor-outdoor exchanges of heat and particulate matter, high resolution multi-layer urban canopy models and building scale Computational Fluid Dynamics (CFDs) are required. Furthermore, for aerosol-climate interactions and microphysics, independent chemistry models coupled to urban climate models are essential. However, this study is about the background perturbations by the synergistic interactions of UHIs and heatwaves (HWs) applied to a regional scale climate assessment, and hence the Weather Research and Forecasting (WRF) model coupled to single layer urban canopy model (UCM) is adequate.

Similarly, heatwaves are one of the most important extreme weather events, because they cause thermal discomfort, heat related health problems, deaths, and higher demands for water and energy resources. Heatwaves are periods of warm weather events that stay for several hours, days, weeks or even seasons. There is no clear nor universal definition of heatwaves. Its meaning varies depending upon the geographical location, environment and the computational methods utilized [5]. An excess of temperature exceeding a certain threshold that stays for two or more days is mostly used. However, such definition cannot be used at a broad regional scale because what is perceived as heatwaves in cooler climates may be taken as a normal hot event in warm climates, i.e., the thresholds (of temperature, humidity and duration) vary regionally [6]. The frequency and duration of heatwaves are also dependent on the state of future climate and human adaptation, hence a warming climate is normally expected to increase the frequency and intensity of extreme weather/climate events, such as heatwaves (e.g., Meehl, [7]).

According to the Inter-governmental Panel on Climate Change’s (IPCC’s) 5^th Assessment Report (AR5), global temperature will increase by more than 2°C if the current rate of emission continues to the end of the 21st century [8]. This will likely result in more intense and frequent heatwaves that would last longer. IPCC also asserted that many lines of evidences show that there is strong, consistent, and almost linear relationship between cumulative $CO_2$ emissions and the projected temperatures at both global [9,10] and regional [11] scales. This indicates that, since most of the emitting populations live in urban regions and modify the local landscape, urban climate is going to be the most vulnerable to the heatwave events. However, if the Paris agreement is successful, the warming will be limited to a range between 1.5°C and 2°C [12].

Nearly every year there were heatwaves in most parts of North America. The human mortality rate due to heatwaves and other extremes of climate such as floods, lightning, hurricanes, tornadoes have been increasing in intensity and frequency [13,14,15,16]. During the periods from 1999 to 2010, over 8081 heat related health problems and deaths had been reported in the US [17]. Most of the affected demographic groups were adults of age ≥65. Almost all the deaths related to heatwaves happened between May and September, with highest numbers reported in July and August. The heat-related deaths occurred with utmost frequency of 81% in urban areas [17]. The number of heat-related deaths in Canadian cities are also increasing [18]. Toddlers and elderlies are the most vulnerable groups. Sensitivity to heat-stress among the elderlies varies from one urban centers to another based on demography, duration of heatwaves, early summer onset, and size of the cities [19].

Conventional assumption dictates that heatwaves are caused by mesoscale and/or synoptic scale atmospheric processes advecting stable, warm and humid air towards susceptible regions [20,21,22,23]. On the other hand, due to urban population growth the need for more buildings, roads and impervious recreational facilities are expected to soar. This in turn requires more land, energy and water, which is expected to exacerbate the urban heat island effects and the consequences of heatwaves. As such, the contributions from local urban centers due to variation of turbulence associated with buildings and surface roughness changes, modifications to the natural thermal and hydrological properties are mostly unaccounted for. But, urban regions are centers of high pressure when under the influence of anticyclones, especially during the night time, when urban surfaces store most of the diurnal solar energy and releases it during the night. This is due to the altered heat capacity and thermal conductivity of urban surfaces, such as asphalt roads, buildings and other impervious surfaces. The phenomenon creates a heat island where the temperature is higher over urban centers and gradually lowers towards the neighbouring rural region in the same microclimate. Regionally, urban heat islands contribute a mean annual surface temperature of around 4°C to the urban climate (e.g., Garuma et al [24]). The impact of such additional urban temperature during the heatwave periods is essential to understand and mitigate the effects of synoptic heat storms. However, there are few studies (e.g., Tan et al, Li et al., Li [25,26,27]) that have linked quantitative contribution of the urban heat islands to the heatwave events based on observations or that were limited to particular and narrow urban areas. The connection between the large scale processes resulting from heatwave events and local scale urban heat islands at intermediate regional scale is poorly studied. Especially, quantitative analyses based on thermal perturbations to the urban energy balance and boundary layer are lacking. The sources of extreme weather events and climate (e.g., heatwaves and coldwaves), and their interrelationships are not well understood [28] in meteorology.

This study, therefore, assesses the synergy between meso-synoptic scale processes, at the regional scale, creating heatwave events and the contributions from local processes using a numerical urban climate model simulations. Numerical urban climate models proved important to study such interactions from microscale to large scale atmospheric processes and the outputs from the Weather Research and Forecasting (WRF) model is used to investigate the interactions between HWs and UHIs employing the responses of urban boundary layer heights, and heat and moisture fluxes. Further analysis was made by linking the heatwave degree days with the urban humidity level and the surface energy partitions. Responses of the surface energy fluxes to a certain threshold of critical temperature for the occurrence of heatwaves are investigated for Eastern North American region. Partitions into pre-heatwave, heatwave and post-heatwave events were performed. These analyses are essential to understand the synergies between HWs and UHIs, and the consequent background perturbations from a larger regional scale perspective.

The aim of this study is to better understand the synergies and interactions of UHIs and HWs. This is useful in identifying the most vulnerable areas to thermal stresses and suggesting ways for the best heat stress reduction strategies. This study provides important information to the future urban planners to assess the impacts of heatwaves and assists prediction and mitigation. Similarly, the synergy information and its subsequent background perturbations can be used to investigate the intensity, duration, and timing of thermal vulnerabilities. Heat related health problems and deaths can be prevented by identifying the most vulnerable areas from a better understanding of processes involved [29].

However, this study is limited to introducing the interaction of heatwaves with urban heat islands over urban surfaces at the regional scale, and not at the detailed urban scale. It deals with investigating the synergies based on the distribution of urban and rural skin surface temperatures, 2 m-humidity, the urban surface energy partitions during, before and after the heatwave event, and responses of urban areas through modifications of regional planetary boundary layer height, surface pressure and potential temperature, and the vertical profiles of potential temperature, humidity and wind. The rest of the paper is organized as follows: section 2 describes the study domain, methodology, initial and boundary conditions and models, section 3 is a detailed discussion of results, and followed by conclusions in section 4.

2. Models, Study Domain, Initial and Boundary Conditions, and Methodology

2.1 Models

In this study, the Advanced Research WRF (ARW) model is applied. It is a widely used community mesoscale model developed by the National Center for Atmospheric Research (NCAR). The WRF modeling system is in the public domain and it is freely available for use by the science community. A detailed description of this model can be found from Skamarock et al [30]. WRF offers multiple options for physics parameterization that can be used in combination with other parameterizations. The options typically range from simple and efficient to sophisticated and more computationally demanding, going from newly developed algorithms to well established schemes such as those used in current operational weather forecasting models. It has wide sets of physical parameterizations available for microphysics, radiation (longwave and shortwave), cumulus and related to the boundary layer: surface layer (SL), planetary boundary layer (PBL), land surface models (LSMs) and urban canopy models of different complexity and details, such as the WRF single layer urban canopy model, WRF-UCM, Multi-layer Building Environment Parameterization (BEP) and Multi-layer Building Environment Model (BEM).

Having all these parameterizations options requires a careful choice for optional configurations. This challenge can be approached in two possible ways, either by performing sensitivity analysis for each possible combination or based on previous studies in science literature. In this study, the best practices from published science literature (e.g. Salamanca et al., Chen et al., Kusaka et al., and Wang et al., [31,32,33,34]) were followed. The most suitable combinations of algorithms were chosen carefully based on the recommendations of previous sensitivity studies that considered urban canopy modeling as a basis for their investigations. As such, the unified Noah [35] land surface model, the single layer urban canopy model, UCM [32], Mellor-Yamada-Janjic (MYJ) TKE scheme [36] for the PBL, Kain-Fritsch (new Eta) scheme [37] for the cumulus physics, Rapid Radiative Transfer Model (RRTMG) [38,39] for shortwave and longwave radiation schemes, and WSM 3-class simple ice scheme for microphysics [40] were selected. The model configuration is summarized in Table 1.

Table 1 Parameterization combinations used in WRF simulation.

The urban canopy model in WRF is similar to the single layer urban canopy model (SLUCM) developed by Kusaka and Kimura, [43]. Urban geometry is represented by two facing building walls separated by infinitely long street canyons. It calculates the heat and moisture fluxes from three different urban surface types given by roofs, roads, and vertical walls. It implements radiation trapping and shadowing effects with an exponential wind profile inside canyons. Heat and moisture fluxes from each of the urban facets are aggregated before transferring back to the WRF model.

The coupling strategy between WRF-UCM and the unified land surface model, Noah is based on the fractional contribution from the natural surfaces (e.g., trees, parks, etc.) and from the urban fraction [32] shown in Figure 1.

Figure 1 The study domain is the entire Eastern North America with the shaded contours showing urban fractions.

The total grid scale fluxes are calculated by the relation

\[ V_{f}=F_{ {nat }} \cdot V_{ {fnat }}+F_{u r b} \cdot V_{ {furb }} \tag{1} \]

where $V_f$ is the land surface heat or moisture flux from the surface, $F_{nat}$ and $F_{urb}$ are natural and urban fractions, respectively.

For the non-urban or rural simulations, where the urban parameterization is switched off, instead a simple treatment of urban land is used. Then a bare ground formulation with modifications to the roughness length, thermal emissivities and conductivities is considered by the unified Noah land surface model. The bulk urban parameterization considered constant parameter values to represent zero-order urban surface effects [44]. Some of the modifications to the PBL to enhance urban surface features are set to increase roughness length representing turbulence generated by surface friction and drag due to buildings. Also, reducing albedo to represent shortwave radiation trapping in urban canyons, altered volumetric heat capacity and conductivity to represent the large heat storage in urban buildings and roads, and reduced green vegetation are applied. However, it does not take into account the effects of urban morphology (e.g., building height, building width, road width), anthropogenic heat, shadowing and radiation trapping effects of buildings. Therefore, the urban surface representation in the Noah land surface model is not an actual urban environment, rather it is closer to bare ground.

2.2 Study Domain, Initial and Boundary Conditions

In this study, Eastern North America (ENA) is chosen because of its relatively higher urban fractions (Figure 1) and its relevance to this sector. It encompasses most of the central, southern and eastern United States, and southeastern Canada. The highest urban fraction in a grid cell reaches up to a maximum of 60% around New England states. Here the urban fraction is the sum of building and road fractions. As shown in Figure 1, most of the highest urban fraction areas are concentrated around the east coast of the United States and Canada. The domain has 184 grid points in the East-West direction and 239 grid points in the South-North direction with a horizontal resolution of 0.22°. Vertically, the standard 62 European Center for Medium-range Weather Forecasting (L62 ECMWF) model levels were used [45].

In this study, the initial and boundary conditions were obtained from the European Center for Medium-range Weather Forecasting (ECMWF) ERA-interim data, with 6h time-step, i.e 0h, 6 h, 12 h, 18 h and 24 h. ERA-Interim is the largest global atmospheric reanalysis produced by ECMWF. Its preparation consisted of input observations interpolated by ECMWF data assimilation system, therefore, ERA-Interim products were used in many meteorological applications [45]. Furthermore, this reanalysis is found to be one of the best performing reanalysis in this North American region (e.g., Sun et al., Beck et al., and Tarek et al., [46,47,48]).

The model is run continuously for the period from January 01, 2000 to December 31, 2002. In this study, we pay a particular attention to the heatwave event of 2000. Most studies that used the WRF model used different spin-up periods (e.g., 24 h by Wang et al., [49] and Lo et al. [50], 12 h by Giannaros et al [51]). Accordingly, in this study, the first 24 h is not used in the analysis to leave time for spin-up period for the model in order to settle to adequate thermodynamic state for removing spurious noise from initial conditions.

2.3 Methodology

Two simulations were performed: one of the simulations consisted of urban surface parameterization and the other did not. For the urban surface parameterization, the single layer urban canopy model, UCM was used. Skin surface temperature (TSK) is selected to quantify UHIs and HWs.

The performance of WRF-UCM in representing the land surface temperature was evaluated using the standard deviation (SD), root mean square deviation (RMSD) and the correlation coefficient (CC) with respect to the Moderate Resolution Spectroradiometer (MODIS) satellite observation data. Taylor diagram [52] was used to evaluate the multiple aspects of the WRF-UCM and WRF models in gauging the relative skill of MODIS. That is, this diagram is used to evaluate whether the models are close enough or further apart in capturing the land surface temperature. The method that produced the land surface temperature closer to the MODIS LST is assumed to perform better.

Once validation is performed and the performance of the models is sufficient, the difference in skin surface temperature (TSK) with a simulation that contains urban surface physics (WRF-UCM) and the simulation without urban surface physics (WRF) is used to represent the surface UHI. Therefore, in this study, UHI is defined as, UHI = TSK_{WRF −UCM} - TSK_WRF. Here the WRF simulated urban surface temperature TSK_{WRF −UCM} is defined as the area average of urban (e.g., roof, wall and road) and rural surface (e.g., vegetation and soil) skin temperatures. The rural skin surface temperature from the WRF simulation (TSK_WRF) uses Noah LSM to simulate the natural land surface fraction (e.g., soil and vegetation).

Heatwave Duration Index (HDI) is used to indicate the extent of heatwaves. HDI is computed as the number of days per time period where in intervals of at least 3 consecutive days, the daily maximum urban temperature exceeds the 90^th percentile threshold [53]. The choice of a relative threshold of temperature is justified since this is a broad regional case study and the absolute threshold temperature for heatwaves varies from place to place [7]. Furthermore, seasonal cycles of UHI, urban canopy humidity, ground heat flux into urban surfaces and urban canopy winds during, pre and post heatwave events are investigated. This provides information concerning which variable controls the heatwave event and its link with the impact of urban surfaces compared to rural. Similarly, the difference between urban and rural surface energy partitions determined from shortwave and longwave radiations, sensible, latent and ground heat fluxes and surface albedos are used to quantify which of these energy processes controls the warming condition during the heatwave event and which factor dissipates it after the event. This could provide clues as to how warming propagates from one season to the other, for example, from summer to autumn. Finally, the impact of the synergies between UHIs and HWs on planetary boundary layer height, surface and 2m-potential temperature is explored in more detail.

3. Results and Discussion

3.1 Validation of the Model Results

The model’s skin surface temperature, $T_{skn}$ outputs were validated against the Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperatures, LSTs. MODIS data is available in 0.05° resolution daily and monthly [54,55] and interpolated to the model’s grid, 0.22° horizontal resolution before comparisons were performed. Because of the improved representation of urban surfaces by WRF-UCM over the WRF only simulation, some improvements were observed due to the expected contribution of the urban canopy processes (Figure 2). It is to be noted that both experiments were carried out using the same initial and boundary conditions and the simulations were performed on the same compute nodes with the same high performance computing (HPC) systems. This is essential to remove the uncertainties associated with the computer resolution and other uncertainties related to inter-model simulations. The statistical performance evaluation using Taylor diagram (Figure 2) was used in this evaluation.

Figure 2 Performance evaluation of WRF UCM (i.e., WRF with urban physics) and WRF without urban physics compared to MODIS observation. The reference satellite data (Ref: MODIS) is shown with the red line. WRF UCM performed better compared to the reference data.

Results compared to MODIS showed that the performance of WRF-UCM is improved over the WRF only simulation (Figure 2). WRF-UCM labelled with a red plus sign in Figure 2 is closer to the reference land surface temperature labelled Ref: MODIS with a red line shown in the diagram. The WRF simulation as shown with a red circle underestimates the land surface temperature compared to the reference LST. Nevertheless, both WRF and WRF-UCM did not deviate much from the MODIS observation. The correlation coefficient is (≈0.98) for both WRF-UCM and WRF simulations. So, the performance of the models is sufficient and can be applied to study the contribution of urban heat islands to the heatwave events.

3.2 Urban Heat Islands and Heatwaves

As explained above, urban heat islands (UHIs) are associated with the material characteristics from which the urban surfaces are built and the consequent energy and water partitions which are different from rural surfaces. On the other hand, a heatwave (HW) is a temperature exceeding a comfortable threshold for a few consecutive days, weeks, or months caused by synoptic scale warming episodes. HWs are mostly consequences of extreme warming periods due to global warming and synoptic events. Understanding the link between UHIs and HWs is essential to accurately predict the vulnerability of urban population to the extremes of temperatures in a warming future climate. Hence, in this study, we used outputs from WRF-UCM and WRF, i.e., a climate model consisting of urban parameterization and without urban parameterizations, to quantify the extent to which both events HW and UHI, are interlinked and the impact the synergy can have on urban population. The differences between outputs from the climate model with the urban canopy and without the urban canopy give an estimation of the extent of the interplay between UHIs and HWs. By how much the urban temperature exceeds the rural temperature during the heatwave events and the consequent responses from the energy partitions and modifications to the urban boundary layer (UBL) give new information whether or not urban surfaces lead to significant retention of the heatwaves over longer periods compared to natural surface (e.g., vegetation and bare ground) and if it favours extension to wider areas, e.g., to suburban and rural neighbourhoods.

Results showed that surface UHI is highest during summer reaching a maximum of 4-5°C (Figure 3) over most of the South Central US. In Autumn, however, its values are reduced to between 3 and 4°C, but the centre of the high UHI values is shifted by about 1° longitudes eastward. Generally, in most parts of the US and some areas of Southeastern Canada, UHI exceeds 2°C. The heatwave event of 2000 happened during summer with a little extension to autumn and hence amplified the UHI during the two seasons (Figure 3). The geographical overlap of the high intensity UHIs and HWs is evidence for a synergistic amplification between HW and UHI. Here the synergy is defined as an increase in urban temperature compared to rural due to the urban heat island during the heatwave period, i.e.,

\[ T_{\uparrow u r b a n, during \ heatwave }=H W+U H I( { synergy }) \tag{2} \]

\[ T_{\downarrow rural,during \ heatwave }=H W- { Cooling }( { non }- { synergy }) \tag{3} \]

Figure 3 Surface urban heat island (Surf-UHI) for each season, (a) winter, (b) spring, (c) summer and (d) autumn.

In an urban environment, the temperature during the heatwave is augmented by UHIs due to the synergistic interaction. This strengthens the intensity of the heatwave over the broader urban regions considering UHI as an additional source of heat. However, the subsequent perturbations to the boundary layer by such additional warming, and energy and water partitions results in a non-linear increase in intensity and duration of the urban temperature. Furthermore, the altered warming boundary could extend the synergistic temperature to suburban and neighbouring rural areas. This has been observed by heatwaves extending to rural areas commonly known as “the warm tongue” [56]. In contrast, the rural temperatures during the heatwave period reduces because of cooling by vegetation and secondary circulations (e.g., lake-breeze circulation and orographic circulations). This is the “non-synergy” between rural surfaces and heatwaves.

It is to be noted that the heatwave happened on the same geographical locations during late summer causing widespread thermal discomfort to people living in the south and southeastern United States.

The maximum heatwave duration index (HDI) overlaps the maximum UHI periods (compare Figure 4 and Figure 3). Most of the locations in the domain exhibited intense and frequent heatwaves in the late summer days, when the heatwave event peaks in August (Figure 5 and Figure 6 (a)). The two other seasons, spring and autumn showed slight warming during exit and entrance heatwave periods respectively. However, the winter season and most days of spring and autumn are too cold to reach heatwave conditions (Figure 5). Noticeably, comparisons of urban and rural skin surface temperatures and the corresponding heatwave events showed that warming is more intense during the heatwave event in urban regions than in rural regions (Figure 6 (a)). This shows that urban surfaces respond more effectively than rural areas to warm periods as indicated clearly by higher surface UHIs in August and September (Figure 7 (a)). This is due to the fact that both urban and rural skin surface temperatures are almost equal at the pre-heatwave and post-heatwave periods, i.e., the differences are higher only during the heatwave event. The differences between urban and rural skin surface temperature increased steadily from June to September (Figure 6 (a) red bars) during the heatwave event. Our simulations show that UHI exacerbates the frequency, intensity and duration of heatwaves. Conversely, heatwaves seem to worsen UHI as exhibited in the same Figure 6 (a), during August and September. This happens because heatwaves do alter the energy balance over urban regions, amplifying the sensible heat fluxes over dry urban surfaces (Figure 8). The differences in emissivities and thermal conductivities of urban and rural surfaces result in uneven energy partitions, which leads to more heat absorbed in urban surfaces and released later during late afternoon or night. Therefore, there is a continuous feedback loop between HWs and UHIs. This supports the idea that during the heatwave event the thermal discomfort is more intense in urban than rural areas. This is discussed further in section 3.3 in terms of the dominant processes involved in the interactions.

Figure 4 Heatwave Duration Index (HDI) calculated based on the mean diurnal temperature exceeding the 90^th percentile for 3 consecutive days.

Figure 5 Number of lat/lon grid points exceeding 90 percentile of the daily mean temperature in the four seasons for the WRF simulation with and without urban canopy model.

Figure 6 Annual cycle histograms of (a) Surf-UHI for urban and rural surfaces, (b)urban canopy humidity, (c) ground heat flux into urban surfaces, and (d) urban canopy wind. Pre-heatwave, heatwave and post-heatwave events are also indicated on the figures.

Figure 7 (a) Urban heat island excess temperature, (b) urban canopy humidity, (c) ground heat flux into urban surfaces, and (d) urban canopy wind for the summer months of June, July, August and September.

Figure 8 The urban and rural energy partitions, and the difference between the two for: (a) net solar radiation over urban surface, R; (b) sensible heat flux, E; (c) latent heat flux, LE and (d) surface albedo, AL.

During the heatwave event, urban humidity is also important. Humidity reduces surface radiative cooling and turbulent mixing by increasing static stability especially during the nights. It is therefore the combination of high temperature and high humidity that cause thermal discomfort. Considering, urban canopy humidity during, pre and post heatwave events of 2000, it peaked during early heatwave event and quickly dropped to lower values at the end of the heatwave event or post-heatwave event (Figure 6 (b)), it decreased steadily from July to October (Figure 7 (c)). This is because urban surfaces dry out more quickly because of the urban impervious surfaces that exposes water to more evaporation during the early heatwave event and then dryness leads to lower humidity at the end of the period. Then, the urban centre completely changes mostly to arid.

Urban canopy ground heat fluxes are indicators of the heatwave events in urban centres because most of the thermal radiations are absorbed at the surface, transferred in urban grounds and released later to flow back into urban centres, which is manifested by the ground heat fluxes. Based on this analysis, urban canopy ground heat fluxes lag by at least one month in responses to the heatwave event, i.e., while the heatwave event peaks during July and August, the ground heat fluxes peak during September (Figure 6 (c) Figure 7 (d)). This is expected provided that urban surfaces have higher heat capacities than rural surfaces, energy is stored into the urban surfaces in summer, especially during the heatwave period and released later in autumn. Therefore, this indicates that considering a warming climate, due to these urban surface characteristics, the heatwave events could be more frequent, intense and have longer duration in large urban regions. This can be inferred from Figure 6 (c)), where the urban canopy ground heat fluxes widened its duration up until January. This shows that in the future, in a warming climate, the summer hot spell condition may extend to fall shortening the cold winter period when considering the heat capacities of the soil and ground and thermal inertia of the built environments.

Wind usually mitigates surface warming through two possible mechanisms, (a) as a natural thermal ventilator, providing cool air from another location to the warm sectors, and (b) providing active weather systems that can disperse the excess of warm air. During dry and windy conditions, obviously the surface temperature lowers appreciably. Dry and cool wind is a natural ventilator that can lower the impact of heatwaves. On the contrary, moist and warm wind may exacerbate the warming conditions. It is because the sensible heat from the warm air increases the local temperature and the moisture may enhance humidity. In contrast, during a heatwave event a reduced wind speed could exacerbate the heat stress. As such, results from this study showed that the wind speed in the urban canopy is reduced on average by 2.5 m/s from January to August (Figure 6 (d) Figure 7 (d)). In addition to the temperature and humidity, wind may also contribute to exacerbating the warming condition in the area. Therefore, it is not surprising that heatwaves could happen when the wind is calm and the PBL is more stable.

In summary, the duration, intensity and frequency of HWs are amplified by UHIs. Urban canopy humidity is higher during the early HWs and dissipates quickly because the urban surfaces gradually run-out of moisture. Urban canopy ground heat fluxes to urban surfaces respond later to the heatwave event and may extend to fall. This would extend the duration of warming in urban surfaces from summer to autumn. Wind speed is lower during the heatwave event over urban surfaces. Therefore, unlike rural surfaces, urban surfaces have optimum characteristics that tend to amplify heatwaves.

3.3 Urban Surface Energy Partitions

The responses of urban areas to heating or cooling events can be determined from the urban surface energy partitions (Figure 9). Knowledge of urban surface energy balance during the extreme weather events could lead to a better understanding of factors mitigating or amplifying urban heat anomalies. These techniques can provide information on how effectively the extremes can be controlled. For example, during heatwaves, should we manage surface albedo, sensible or latent heat fluxes via vegetation or drainage from urban surfaces to achieve optimal urban cooling? To answer this question, first we need to interpret the urban energy balance model and then identify how it deviates from that in rural areas. Later, the knowledge can be used to relate it to its responses during a heatwave event and the information can be used to design the best cooling strategies in urban environments.

Figure 9 Illustration of the energy balance fluxes and boundary layer heights during heatwave and non-heatwave periods for building-vegetation-soil-air model in a grid cell. The energy balance terms E^∗,R, L, E, LE, E_f , E_c, E_a, are the net shortwave radiation, the incoming solar radiation, longwave radiation, the turbulent sensible and latent, anthropogenic, the ground conduction and advection heat fluxes. The vertical height, h₁ is at the urban canopy layer height (UCLH), h₂ is at the planetary boundary layer height when there is no heatwave (PBLH_{no−heatwave}) and h is when there is heatwave (PBLH_heatwave).

The energy balance equation (e.g., Oke [57]) for an urban surface is

\[ Q \ *=R+L=E+L E+E_{N} \tag{4} \]

where $Q_*$ is the net radiation on urban surface, R is the net solar radiation, L is the net infrared radiation, E is the sensible heat flux and $E_N=E_f+E_c+E_a$ is the net of other heat fluxes: anthropogenic, conductive (ground) and advective heat fluxes. The net solar radiation in urban microclimate is affected by urban geometry, i.e., the position, orientation and shapes of buildings. The net solar radiation at the surface is

\[ R=\left(1-\alpha^{-}\right)(D S+R R+D D) \tag{5} \]

where $\alpha^{-}$ is the surface albedo, DS is the direct solar radiation received by the surfaces, RR is the solar radiation reflected by the surface j (e.g., the left wall) and received by surface k (e.g., right wall) and DD is the diffuse solar radiation received by the surface. The diffuse and reflected solar radiations components are functions of the sky view factor and the relative view of one surface type (e.g., wall) with respect to the other (e.g., road).

In an urban settings, due to attenuation by urban particulate matter, the radiation reaching the urban surfaces are lower than that over rural surfaces. This is the case when there are significant amounts of aerosols in the atmosphere. On the other hand, the net long wave radiation is higher in cities due to blockage by urban sky view factor, i.e., more long wave radiation is trapped in urban canopies than released into the free atmosphere. Therefore, the net radiation over urban surfaces is almost in balance. However, this could vary depending upon the microclimate of the urban area, i.e, the surface energy balance varies under clear sky and overcast conditions in addition to the distributions of air pollutants. On the other hand, perturbations due to extreme weather events, such as heatwaves, could lead to changes to the extent of disturbing the energy partitions, boundary layer depth and vertical mixing.

During the heatwave event, the difference in net shortwave radiation between urban and rural surfaces is generally low (from -3 W.m⁻² to 1 W.m ⁻², Figure 9 (a)). The shortwave radiation is lower in all the months except during the heatwave months (e.g., July and August) (Figure 9 (a)). The source of the imbalance is the result of multiple reflections and diffused radiation (eqn. 5) by aerosols downwards to urban surfaces. A very small change in net shortwave radiation could lead to more urban surface heating.

Considering longwave radiation over urban and rural surfaces, the difference is as expected higher over urban surfaces than rural during the heatwave period (Figure 9 (b)). However, its intensity remains higher for a few months after the event. In July and August, both net surface shortwave radiation and longwave radiation over urban surfaces is higher than over rural surfaces. This clearly indicates that urban surfaces worsen the heatwave compared to the rural surfaces. Surface albedo, $\alpha^{-}$in eqn. 5, is quite different in urban and rural environments and significantly affects the net energy balance over urban surfaces. As shown in Figure 9 (e), the surface albedo is lower during summer season, especially during the heatwave event in July and August. So, much of the incoming solar radiation is absorbed by urban surfaces and exacerbated the warming due to the warm air masses. Lower surface reflectivity followed by high absorptivity over urban surfaces are characteristics that are the causes for urban heat islands and these amplified the intensity, duration and coverage of the heatwave event. Furthermore, as urban surfaces are highly impervious and have less vegetation, there is less water left for evaporative cooling as exhibited by lower latent heat fluxes, LEs (Figure 9 (d)) during the heatwave event. The latent heat flux is much lower during the peak of heatwave events in July, August and September resulting in lower evaporative cooling in urban areas, which further increases the duration and intensity of heatwaves.

On the other hand, since the energy conversion to other forms (e.g., latent heat) is lower, most of the absorbed energy by urban surfaces, Q* is converted to sensible heat flux, E (Figure 9 (c)). During the heatwave event, in summer both the net solar radiation (the vector sum of shortwave and longwave radiation) and sensible heat fluxes are higher. This shows that most of the heatwave events are followed by higher net radiation augmented by turbulent surface characteristics. In summer, the energy difference between sensible heat and net solar radiation is absorbed by the surface and partly converted to latent heat. The higher sensible heat flux and net solar radiation during the heatwave event over urban surfaces signify that the urban heat island worsens the heatwave event over urban surfaces. This suggests that both urban albedo and evaporative cooling management (e.g., cool roofs and vegetation) are viable means of reducing future heatwave events.

In conclusion, during the heatwave event the urban surfaces favor more urban heating than cooling, which is exhibited by higher solar radiation, sensible heat fluxes, ground heat fluxes to urban surfaces, and lower surface albedo and latent heat fluxes. Over urban areas, therefore, the energy balance tends to amplify heatwaves. This is not the case for rural surfaces, where vegetation and soil moisture favor cooling the heatwaves compared to the urban surfaces where the heatwave is amplified.

3.4 Boundary Layer, Surface Pressure and Potential Temperature

The lowest part of the troposphere influenced directly by the Earth’s surface through exchanges of heat, momentum and moisture is the planetary boundary layer (PBL) [58,59]. PBL has an effect on the surface energy budgets because it modifies cloud properties, coverage and the mixing heights of pollutants. The exchanges of moisture, heat, and momentum within the PBL occur through mixing associated with turbulent eddies. The evolution of the lower-tropospheric thermodynamic and kinematic structures is governed by these eddies. However, most mesoscale models cannot explicitly represent the spatio-temporal scales in which the eddies operate. The effects of the eddies are expressed in the models through the use of PBL parameterization schemes. However, much less is known about the vertical extent of the urban-induced modification [60]. The Mellor Yamada Janjic (MYJ) Scheme used in this study determines the PBL height (PBLH) using the TKE profile. It estimates the PBL top at the height where the TKE decreases to a prescribed low value (Janji´c, [42] as reported by Hu et al. [61]). It is necessary to determine PBL characteristics during heatwave events to investigate the nonlinear relationship between urban surface characteristics and its influence on energy and water partitions in the boundary layer and compare it with the rural boundary layer.

During heatwave events, the depth of PBLs are heightened (Figure 9), especially in urban areas (this is inferred from the difference in PBL height over urban and rural areas, which is a positive number), where the influence of the heat islands is higher and more of the surface energy fluxes are from sensible rather than latent heat.

The mean planetary boundary layer height (PBLH), H, is the result of mechanical and thermal surface characteristics, and during normal atmospheric conditions, it is

\[ H=\left(H_{m}+H_{t}\right)_{U, R} \tag{6} \]

where $H_m$ and $H_t$ are the mechanical and thermal contributions for urban and rural regions represented by subscripts $U,R$. During this normal event, the urban boundary layer height is greater than rural, i.e., $\left(H_{m}+H_{t}\right)_{U}>\left(H_{m}+H_{t}\right)_{R}$. This is because of the mechanically induced urban surface lengths that are larger than rural surfaces, and the thermally induced perturbations by the urban heat island.

During extreme climate events of heatwaves and cold waves, an additional thermal perturbation $H^{\prime}$ leads to

\[ H=\left(H_{m}+H_{t}+H^{\prime}\right)_{U, R} \tag{7} \]

The perturbed factor, $H^{\prime}$ is negative for coldwaves and positive for heatwaves, because the cold period lowers the PBLH, while the warm period heightens it. Because of the urban thermal characteristics, the urban response to thermal heating is to impart an additional PBLH, i.e., $H^{\prime}{}_U>H^{\prime}{}_R$ .

In the urban air, the shift from latent to sensible heat flux is associated with deepening PBL, $H^{\prime}{}_U>H^{\prime}{}_R$ [62]. Our results confirm this fact with the PBL being higher during the heatwave months, July and August (Figure 10). The peak of the PBL height difference between urban and rural ($H_{U-R}$) reaches up to 90 m from the ground, which appears mostly during the months of July and August (Figure 10 (i)), i.e., PBL responds to the urban surface heating with equal proportion but later during the heatwave period. This is because of the response differences in heat capacities of the urban air and morphology compared to rural. The unequal response from heating by heatwaves and urban heat islands creates a temperature gradient optimum for the occurrence of free convection. Convection is an important means of diffusing heat and pollutants into large atmospheric volumes, and the depth of the mixed layer sets an upper limit to the vertical dimensions of this volume. Consequently, the increase in the vertical PBL pushes the limit in the horizontal direction. Therefore, the effect of deep PBL height is that it increases the surface areas vulnerable to the synergistic heating by the urban heat islands and heatwaves. This can extend the tongue of the UHI (known as “urban heat plume”) towards neighbouring rural areas. In a warming climate, a dry surface heats up more quickly than a moist surface, which causes more longwave radiation flux, reducing the available energy at the surface [63]. Hence, urban surfaces are associated with reductions in the surface turbulent energy exchanges between the surface and the planetary boundary layer. The energy deficits would create a laminar boundary layer immediately above the surface. So, the adjacent layers of the fluid remain distinct and do not intermix. This smooth flow would in turn reduce the mixing in the PBL and intensify the near surface heating during the heatwave event as shown by the peak of the potential temperature in August and September (Figure 11) and the higher near surface temperature difference between urban and rural areas (Figure 12 (c)). The thickness of the laminar boundary layer reduces gradually until turbulent flow takes over at a critical height above the surface of the laminar flow layer. In the thin laminar surface near the ground there is no convection, therefore, all non-radiative transfer is by molecular diffusion, which is inefficient for mixing. The heat trapped in the thin laminar flow near the surface stays there for a long time. In contrast, immediately above the laminar surface into the turbulent eddies above the laminar flow, leading to enlargement of the surface area with thermal transport.

Figure 10 Monthly mean PBL height for urban, H_U (column 1), rural, H_R (column 2) and the difference, H_U−R (column 3) during the heatwave period (June, row 1; July, row 2; August row 3; and September row 4). The first label bar to the right is for the urban and rural PBLH and the second label bar to the left is for the difference between urban and rural PBLH values.

Figure 11 Monthly mean potential temperature during the heatwave period (June, row 1; July, row 2; August row 3; and September row 4) for urban (column 1), rural (column 2) and the difference (urban - rural) (column 3). The first label bar is for the urban and rural and the second label bar is for the difference between urban and rural values.

Figure 12 Monthly mean vertical profiles of temperature (first row) for, (a) urban, T_U, (b) rural, T_R and (c) the difference, T_U−R, humidity (second row): (d) urban, q_U, (e) rural, q_R and (f) rural, q_U−R, u-wind (third row): (g) urban, U_U, (h) rural, U_R, and the difference, U_U−R(i) and v-wind (fourth row): (j) urban, V_U, (k) rural, V_R, and (l) the difference, V_U−R.

Multiple elevated inversions over the cities (warm air from above and cold air from below) extend the top of the boundary layer (Figure 10). Therefore, the response of the PBL is such that it prolongs the heatwave period (because of the heat trapped in the laminar boundary layer), extends the areas covered by the UHIs (because of the deep turbulent layer above the laminar layer) and intensify the near surface heating (because of the solar radiation trapped in the deep mixing layer). Considering surface pressure, it is lower when the PBL height is higher earlier during the heatwave (June and July) than later (August and September) (Figure 13). Specifically, the higher pressure in the urban simulations (Figure 13 (f)) is replaced by lower pressure in August (Figure 13 (i)) over Eastern Canada and Southeastern US. This implies the heavy and colder air during July is replaced by warm and convective air during August. This shows that urban regions can alter the regional climate system. However, the magnitude of the pressure change is very low and hence less correlation is expected. Nevertheless, early surface heating increases the surface pressure and the gradient disappears towards late summer, which signifies, there is small surface pressure change associated with the heatwave and the consequent urban heat islands.

Figure 13 Monthly means surface pressure during the heatwave periods (June, row 1; July, row 2; August row 3; and September row 4) for urban (column 1), rural (column 2) and the difference (urban - rural) (column 3). The first label bar is for the urban and rural and the second label bar is for the difference between urban and rural values.

3.5 Sensitivity to Vertical Mixing

Vertical mixing of temperature, water vapor and wind during heatwaves can explain the level of perturbations to the PBL dynamics by the land surface heating. The synergistic heating of the land surface by UHIs and HWs further extend the perturbations to a higher level over urban land surfaces. Results show that the profile of air temperature in the lowest 5 km of the atmospheric profile is higher during the heatwave event (months of July and August) for both urban and rural areas (Figure 12 (a) and Figure 12 (b)) than the other months. However, near the ground, up to 1.0-km from the land surface, the urban temperature exceeds the rural ($T_{U-R}$) by up to 0.12°C (Figure 12 (c)) indicative of the synergy of UHI and HW. Due to the deepening of the PBL during the heatwave event (refer to section 3.4 and Figure 10), the thermal convection further protrudes from the mixing atmospheric layer to the turbulent atmospheric layer - due to the PBL thermal inertia. However, above 1 km, the difference in urban and rural temperature, T_U−R becomes uniform during the heatwave period (refer to the red line for the month of August, Figure 12 (c)). This is when the urban influence is reduced at a distance further from the urban canopy layer (UCL). This additional near the ground temperature, T_U−R could create a thermal inversion over urban surfaces during the late night. Therefore, the interaction of HWs and UHIs favors the warming of the lowest atmosphere exacerbating thermal discomfort to the urban population during the day and the formation of thermal inversion during late night or early mornings. Subsequently, this would create lower thermal mixing that can have serious weather and climate consequences, for example, pollutant stagnation.

Along with temperature, humidity is an important factor, especially during the heatwave event (Figure 12 (d) and Figure 12 (e) red lines). High temperature augmented by water vapor causes thermal discomfort, heat-related illnesses or death. Typically, urban canopy humidity (UCH) is higher during early heatwave period and tends to dissipate afterwards (section 3.2, Figure 2 and Figure 7) because less urban surface water is available for evaporation in the later phase of HW and UHI (note that the water vapor mixing ratio, $q_{_{U-R}}$ is higher during early heatwave events, June and July and decreases thereafter, Figure 12 (f)).

The other important parameter for the synergistic UHI and HW interaction is wind. As explained in section 3.2, wind (both urban canopy wind and the vertical profiles of wind) were extremely low during the heatwave event (Figure 12 (c)). The surface temperature peaks during July and August (Figure 12) but U-wind speed is extremely low during these months. The low ventilation and little mechanical turbulence due to weak air circulation create a relatively stagnant air over the area during the heatwave event. Such low wind speeds provide little ventilation that permits the heatwaves to stay longer. Considering the vertical profiles of wind, both the U-wind (Figure 12 (g) to (i)) and the V-wind are lower for the urban areas, i.e, $U_{U-R}\leq 0$ and $V_{U-R}\leq 0$ during the heatwave event (July, August and September). This shows that urban areas tend to block or constrict the free air circulation and this would create a calm weather condition that would be favourable for further warming. In contrast, rural areas favor cooling by evapotranspiration and the relatively less obstacle to the free air circulation.

Furthermore, understanding the level of perturbations by the synergistic interactions of UHIs and HWs can give us background information useful for mitigating heatwaves. PBL parameterization schemes are based on the decomposition of variables of the equations of motion into mean and perturbed components. The background atmospheric state is represented by the time-averaged mean components of the equations. On the other hand, the perturbation components represent deviations, or turbulent fluctuations, from the background mean state. The turbulent kinetic energy (TKE) is used to quantify the perturbation component of momentum in the Mellor and Yamada (MYJ) scheme [36]. It provides the extent of turbulence in the PBL with respect to vertical wind shear, buoyancy, turbulent transport, and dampening driven by molecular viscosity (Holton et al. [64], as reported in Cohen et al. [65]). Because the MYJ PBL scheme uses TKE to determine the PBL height, TKE (Figure 14) shows similar profiles to the PBL height. Near the ground, upto 1.5 km, the TKE_U−R is around 1.5 × 10⁻³ m⁻².s⁻², during the months of August and September (Figure 14 (c)). Analogously, the length scale from PBL (Figure 14) is highest during the month of September (Figure 14 (f)). This implies that the total atmospheric kinetic energy over urban areas extends to late summer because of the urban surface characteristics, i.e., urban surfaces can absorb, retain and extend the impact of the synergies (UHI+HW) to an extended period of time.

Figure 14 Seasonal cycle of Turbulent Kinetic Energy (TKE) for (a) urban, TKE_U (b) rural, TKE_R and (c) the difference (TKE_U−R), and the length scale from PBL for (d) urban, l_U (e) rural, l_R and (f) the difference (l_U−R).

4. Conclusion

In this study, validation of the land surface temperature (LST) from the simulations of the WRF model, with and without urban canopy parameterization, was performed against high resolution MODIS observation data. Using the validated results, the synergistic interactions between HWs and UHIs, and the subsequent background perturbations to urban surface energy balance and planetary boundary layer height were evaluated. Results showed that during the heatwave event urban surfaces enhanced the UHI which emanates from higher solar radiation due to lower surface albedo, higher sensible heat flux, ground heat fluxes to urban surfaces and lower latent heat fluxes than the surrounding rural areas. UHI contributes a maximum seasonal mean land surface temperature of around 4°C during summer (when the heatwave happened). The intensity of both UHI and HW overlap geographically which is evidence for a synergistic amplification of the urban surface heating. The synergy, in this study, is defined as an increase in urban temperature (due to UHIs) during heatwaves compared to rural areas, i.e., in an urban environment, the temperature during HWs is augmented by UHIs. Heatwaves do alter the energy balance over urban regions favouring more sensible heat over latent heat flux, and this has exacerbated the heatwaves over urban regions than rural. This is partly attributed to the differences in thermal emissivities and conductivities between urban surfaces (higher impervious land fractions) and rural surfaces (higher pervious land fractions and vegetation). This is not the case for rural surfaces, where vegetation and soil moisture favor cooling the heatwaves compared to the urban surfaces where the heatwaves are amplified.

In urban air, the fractional shift from latent heat to sensible heat flux is associated with deepening of the planetary boundary layer so that when the height of the PBL is higher during the onset of the heatwave event or during the pre-heatwave period, the duration and extent of the heatwaves increases. Humidity also showed substantial shifts during, pre- and post heatwave period. It is to be noted that humidity is an important meteorological parameter during heatwaves because humidity reduces surface cooling by reducing turbulent mixing in the planetary boundary layer. Our results indicated that urban humidity decreased steadily from pre-heatwave to post-heatwave implying urban surfaces dry out quicker than rural areas because of the more impervious urban fractions that exposes water to more evaporation during early heatwave event and then dryness leads to lower humidity at the end of the heatwave period, at which point the urban center completely changes to aridity. The other important parameter that influences the surface heating through thermal mixing during heatwaves is wind.

Wind is a natural cooling agent, i.e., it transports cool air from a distant location to the warming area, simultaneously providing a pressure force that would disperse the warm air. Results from this study showed that wind in the urban canopy was lowered by 2.5 m/s from pre-heatwave to heatwave event. This shows that the air dynamics is very limited resulting in no substantial impact on heatwaves. Comparing urban and rural regions, the wind speed in the urban canopy dropped by 0.15 m/s for urban areas compared to rural. Low wind speed tends to extend the heatwave period, reducing the flow of cool air towards urban areas and less pressure gradient to disperse local warm air packets. Therefore, in conjunction with temperature and humidity, wind may also have its own contributions in exacerbating the warming condition in urban areas. As such, it is not surprising that heatwaves worsen when the wind is calm and the PBL is more stable.

Determining the PBL characteristics during heatwaves is important to investigate the nonlinear relationship between urban surface characteristics and its influence on energy and water partitions in the urban boundary layer (UBL) and compare it with the rural boundary layer (RBL). The responses of the PBL in urban areas is in such a way to prolong the heatwave period, extending the areas covered by the UHIs (because of the heat trapped in the laminar boundary layer) and intensifying the near surface heating (because of the trapped solar radiation near the ground). The shift from latent heat to sensible heat flux in urban air is associated with deepening the PBL. Our results confirm this such that the peak of the PBL height difference between urban and rural ($H_u-H_r$) reaches up to 90 m from the ground up, which appeared during the heatwave months of July and August. As such, the increase in PBLH by the synergistic UHIs and HWs tends to widen the surface areas vulnerable to the combined heatings. This widening could stretch out to neighbouring suburban and rural areas.

Vertical mixing of temperature, humidity and wind show that due to the deepening PBL during the heatwave event, the thermal convection protrudes further to higher altitudes over urban regions than rural, i.e., up to 90 m from the surface. Above this altitude, the temperature drops quickly during the heatwave event in summer more than the other seasons. Nevertheless, knowledge of the level of atmospheric perturbations by the UHI+HW synergies is important to get information useful for mitigating heatwaves. It is to be noted that the background atmospheric state is the sum of the time-averaged mean state of the equations and the perturbed components. The deviations, or turbulent fluctuations from the background mean state is represented by the turbulent kinetic energy (TKE). In this study, TKE_U−R showed similar profiles to the PBLH. The difference in urban and rural turbulent kinetic energy, TKE reaches around $1.510^{-3} \mathrm{~m}^{-2} \cdot \mathrm{s}^{-2}$ during and post-heatwave period implying the total TKE over urban areas extend to late summer where the urban surface characteristics retain the UHI and HW synergies to an extended period of time. This study pointed out the best urban heat-stress management strategies from the perspectives of urban energy balance and vertical mixing in the planetary boundary layer. Thus, increasing surface albedo and latent heat fluxes, lowering sensible heat fluxes through urban vegetation and increasing the flow of cool air towards warm urban centres (the urban ventilation corridors) are the recommended methods to reduce heat related illnesses and mortalities in the future warming climate.

Acknowledgments

All the simulations and data analyses considered in this study were performed on the supercomputer managed by Compute Canada and Calcul Quebec. I would like to forward special thanks to the Cedar heterogeneous cluster system management team because most of the simulations were performed on their cluster computer.

Author Contributions

The author did all the research work of this study.

Competing Interests

The author has declared that no competing interests exist.