Urban Heat Island Dynamics in an Urban–Rural Domain with Variable Porosity: Numerical Methodology and Simulation

: Heat transfer and ﬂuid dynamics modeling in porous media is a widely explored topic in physics and applied mathematics, and it involves advanced numerical methods to address its non-linear nature. One interesting application has been the urban-heat-island (UHI) numerical simulation. The UHI is a negative consequence of the increasing urbanization in cities, which is deﬁned as the presence of warm temperatures inside the urban canopy in contrast to the colder surroundings. Furthermore, an interesting phenomena occurs within a UHI context when the city transitions from a heat island to a cold island, matching the increases and decreases of solar radiation over the span of a day, as well as the decrease in the UHI intensity as a result of wind action. The numerical study in this paper had, as its main goal, to reproduce this phenomenon. Therefore, the key elements proposed in this work were the following: A 2D horizontal urban–rural domain that had a variable porosity with a Gaussian distribution centered in the city center. A non-stationary Darcy–Forchheimer–Brinkman model to simulate the ﬂow in porous media, combined with an air–soil heat transport model linked by a balancing equation for the surface energy that includes the evapotranspiration of plants. In regards to the numerical resolution of the model, a classical numerical methodology based on the ﬁnite elements of Lagrange P 1 type combined with explicit and implicit time-marching schemes have been effective for high-quality numerical simulations. Several numerical tests were performed on a domain inspired by the metropolitan region of Guadalajara (Mexico), in which not only the temperature inversion was reproduced but also the simulation of the UHI transition by strong wind gusts.


Introduction
The urban heat island (UHI) is a well-known phenomenon that outlines the discrepancies between the temperatures of an urban region and the surrounding rural region. The UHI itself has a different behavior between day and night, being more remarkable at night and during the summer season [1][2][3], in contrast with the complex distribution in the daytime with warmer areas depending on soil use (industrial, parks, commerce, residential, and others) [4]. Moreover, the exposure of city inhabitants to these warmer temperatures generate thermal stress, defined as discomfort and feeling upset with the consequences of the warmer temperatures on their daily activities, such as working and sleeping, and in worst-case scenarios, these have even resulted in death due to heat waves [5,6]. In order to mitigate the thermal stress, an extensive use of cooling systems is the only alternative, leading to increased energy consumption in buildings by a median rate of 19% [7]. Furthermore, the UHI facilitates the presence of higher pollution levels in cities. This was incompressible turbulent flow (Forchheimer correction), along with the fluid that tended to adhere to the walls of the porous media (Brinkman correction). We could infer the reasons for this contrast, as an entire urban-rural domain is unlikely to be supported by the experimental results of the above models. Regarding the heat transfer, [12,16,20] focused on the air temperature (fluid phase) forced by a known heat flux, depending on the building height and volume, enabling them to eliminate the equation for the solid phase, which complemented the heat-transfer model on porous media [21]. Therefore, our work identified remarkable contrasting elements, with respect to [12,16,17,19,20], which are mentioned here briefly and are detailed throughout the paper. Our urban-rural domain is a 2D horizontal region with a variable porosity defined as a Gaussian distribution to focus on a human scale (a few meters above the soil surface), overlooking the turbulent vertical transport. Furthermore, similar to [20], we apply a non-stationary CFD model of a Darcy-Forchheimer-Brinkman type and coupled with a heat-transfer model that includes not only the fluid-solid phase equations found in the porous media literature [21] but also an equation for the energy balance on the surface that has been used in remote sensing research [22]. Therefore, our model has solar radiation as the energy source, instead of heat fluxes from data or buildings. Lastly, we show the complete mathematical methodology and the numerical schemes used to solve the two models and carry out several numerical experiments. All of this is considered with the aim to reproduce the phenomenon of temperature inversion and the wind transport of warm air. This paper has the following structure: In Section 2.1, the elements of a typical domain are described, together with a brief justification of the model. In Section 2.2, we present the momentum equation as a Darcy-Forchheimer-Brinkman equation but re-written as an equation depending on the local pore-scale in the representative elementary volume (REV), as well as stating our motivation and its advantages. Next, the heat-transfer model is presented in Section 2.3, in which we highlight the inclusion of an equation for the energy equilibrium on the surface that considers the total solar radiation, the heat fluxes of sensitive and latent types, and the soil. In Section 2.4, the finite element method and timediscretization schemes for solving both models are shown. While the numerical schemes are classic, modifications due to variable porosity were necessary, and boundary conditions had to be completely described for the purpose of presenting a complete numerical methodology, which has been absent in much of the literature, such as in [12,16,19,20]. In Section 3, the model domain is set as the metropolitan region of Guadalajara (MRGDL) in Mexico, and it is characterized not only by a variable porosity but also by other parameters, such as density, specific heat, albedo, and emissivity, among others. Numerical experiments are carried out and show the UHI evolution across a 24 h period and the transport of warm air mass due to wind gusts. Finally, in Section 4, we explain how the model achieved the expected results by simulating the urban temperature inversion and the wind action over a UHI.

The Computer Fluid Dynamics and Heat Transfer Models on Porous Media
Currently, heat transfer and fluid dynamics through porous media are phenomena studied by physicists and applied mathematicians, where the numerical simulation of PDE models is an essential tool. It is our interest to study the corrections to Darcy's Law, which consist of added terms to address fluid behaviors observed in laboratory experiments, such as those for high-velocity flow and resistance on the porous wall (viscosity). These new terms or corrections are known as Forchheimer and Brinkman terms, respectively. In regards to modeling the heat transfer, two phases are considered: fluid phase and solid phase, both formulated as typical non-stationary transport problems with advection, diffusion, and averages in vertical forcing terms. Furthermore, fluid and solid phases are linked by a heat-interchange surface that plays the role of a temperature regulator. Therefore, an indispensable parameter that is present in both models is the porosity, which is dimensionless with values between 0 and 1, indicating the proportion of pores in the REV. As mentioned above, a UHI is characterized by the contrast between temperatures of urban and rural media; therefore, we consider a regular domain Ω ⊂ R 2 with boundary Γ divided into different segments: wind inlet Γ in , wind outlet Γ out , and hills inside the domain Γ w .

A Non-Stationary Darcy-Brinkman-Forchheimer Model for an Urban-Rural Porous Media
Let Ω ⊂ R 2 be a porous media with porosity = (x) and permeability K = K( ). Let (0, T) be the time interval of the simulation. We had to find the average velocity on the reference elementary volume (REV) u = (u 1 , u 2 ), such that it satisfied a CFD Darcy-Forchheimer-Brinkman model, formulated as follows [21]: for index k = 1, 2 and where P is the pressure, ρ is the air density, µ is the dynamical viscosity, C F is the Forchheimer coefficient, and u also satisfies the continuity equation Let the local fluid velocity on the REV be defined by u * = u , so the above equation can be expressed in terms of u * , as follows: The left hand side of the last equation is well known, and this method for defining u * allow us to send the porosity = (x) to the forcing term on the right hand side of the equation, which is more advantageous when seeking its numerical solution. Accordingly, the continuity equation needs an adjustment, as a consequence of the variability of the porosity. With respect to the initial and boundary conditions, we assume a known initial state for the velocity field. Therefore, for a given u 0 (·), Meanwhile, the boundary conditions are formulated based on the boundary segments, as follows: where n = (n 1 , n 2 ) is the unit normal vector at each point of Γ and exterior to the domain Ω. The boundary conditions indicate that a constant inlet wind is imposed on the inlet boundary Γ in (6a), and a no-slip condition is used on walls Γ w (6b). Furthermore, boundary condition (6c) indicates that the wind direction on Γ out is not affected by the pressure. Therefore, pressure on the inlet boundary Γ in and walls Γ w is applied in a tangential direction under condition (6d).

Model for Heat Transfer in Urban-Rural Domain
We are interested in simulating the temperature dynamics for our urban-rural domain. In order to achieve this objective, three temperatures are considered: air, surface, and soil, with solar radiation as their energy source. Next, air temperature is important for quantifying the UHI intensity, and it is strongly influenced by the wind flow field. The surface acts as a screen for the solar radiation, reflecting a fraction of it and absorbing its complement. As a consequence, the surface temperature can be determined. Meanwhile, soil temperature depended only on the energy interchange with the surface. Lastly, despite being defined on the 2D horizontal domain Ω, the mentioned temperatures are related by heat transfer at a much smaller vertical component.
Let θ a , θ 0 , and θ s denote the air, surface, and soil temperatures, respectively, which satisfy the following PDE system: where we observe the porosity , the global total radiation R s , and those related to air, such as µ a , γ a , d a , and σ a , as the thermal diffusivity, convective interchange coefficient, air layer thickness, and radiation transfer coefficient, respectively. Parameters related to the soil are µ s , γ s , and d s which correlate to the thermal diffusivity, convective interchange coefficient, and soil layer thickness, respectively. Furthermore, those related to the air-surface-soil link are the albedo a, the sky emissivity sky , the soil emissivity 0 , soil interchange resistance r sh , and air interchange resistance r ah . The complete list of the model parameters is found in Table A1.
For this system, we identified the following elements: Equation (7b) represents the energy balance on the surface, commonly used to estimate surface temperatures or other surface variables in remote sensing [22,23]. This equation formulates the equilibrium between total radiation , and latent heat flux λLE. The latter is obtained from the definition of the Bowen ratio, β = H λLE , and the energy balance equation on the soil surface: Then, substituting the value H = βλLE, we obtained a method for estimating the latent heat flux with the following relation [24,25]: as long as a known value of β is available.

Remark 1.
Additional physical definitions are required to explain Equation (7b). The total radiation is the difference among the downward short-wave radiation or solar radiation (1 − a)R s , the longwave radiation from the sky sky σ B θ 4 a , and the upward long-wave radiation emitted by the soil 0 σ B θ 4 0 . The sensitive heat flux H and latent heat flux λLE are the mechanisms by which the surface transfers much of the absorbed solar radiation to the air. The magnitude of H depends on the temperature difference and the resistance r ah while λLE depends on the Bowen radius, which is related to evapotranspiration, predominantly by plants. Therefore, H and λLE play the role of thermal regulators for the air layer temperature.
The Equations (7a,c) are formulated as a typical transport problem with diffusion, advection, and vertical averaged forcing, by d a and d s , respectively. Furthermore, the location on the domain of the three temperatures and its relation with the convective resistances r ah and r sh is shown and briefly explained in Figure 1. Figure 1. Heat-transfer model: the urban-rural domain Ω (bottom) is a sub-set of R 2 with an urban area, represented as a rectangle (gray); a rural area (green); the soil below (brown); inlet-outlet boundary segments; and inlet wind direction (blue arrow). Based on the horizontal section of the domain (magenta sector), the vertical distribution of the three temperatures θ a , θ 0 , and θ s is shown (top left) in congruence with the vertical interchange thermal resistances r ah and r sh for air and soil, respectively (top right). Both layer thicknesses, d a for air and d s for soil, are much smaller when compared with the horizontal scale. Finally, the resistance diagram (top right) is formulated based on Figure 6 in [23].
The system is complemented with the initial conditions for air and soil temperatures θ a (·, 0) = θ 0 a (·) and θ s (·, 0) = θ 0 s (·), where θ 0 a (·) and θ 0 s (·) are given functions. Regarding boundary conditions, it is assumed that heat escapes from the domain only by wind action in the case of the air temperature, so diffusion must be at zero on all borders. Whereas for the soil, it is assumed that soil temperature is equal on both sides of the boundary, implying zero diffusion. Therefore, the following homogeneous Neumann boundary conditions for both air and soil equations are imposed: where n is the unit normal vector at each point of Γ and exterior to the domain Ω.

Remark 2.
Conditions (10a,b) have the advantage of simplifying the variational formulation of the model, and both conditions are a special case of the heat interchange condition n · µ∇θ = −φ(θ re f − θ), where φ > 0 and θ of the model are unknown and θ re f is a given reference temperature.

Numerical Solution
Notice that the air density is constant, so the systems (3) and (7) are uncoupled, and they can be addressed independently. Moreover, the Darcy-Forchheimer-Brinkman Equation (1) has a viscous term that contributes to stability; meanwhile, the equations of the energy transfer model had the classic form of a transport equation. This indicates that a standard finite element methodology can be applied. Therefore, let us consider the Sobolev space V = H 1 (Ω) as our space of test functions, along with the L 2 inner product ·, · to obtain the variational formulation, and the finite element functions of Lagrange P 1 type to achieve a stable numerical solutions [26][27][28].

Remark 3.
Note the mathematical analysis of the Darcy-Forchheimer-Brinkman and Darcy-Forchheimer equations in order to define the appropriate Sobolev spaces and to derive suitable variational problems has been studied in [29,30], thereby supporting the use of a finite element approach. However, due to the scope of this work, it is sufficient to use a variational approach based on the space H 1 (Ω).

An Explicit Scheme for the Darcy-Forchheimer-Brinkman Equation: The Chorin Method
The Chorin method is now considered a classic method for solving the Navier-Stokes equations numerically [31], which is explicit in time and subsequently splits the model into three equations. The method has the advantage of being straightforward to formulate, at the expense of losing accuracy on the boundary conditions [26]. The method begins by formulating a semi-discrete problem in time. Let T denote the final simulation time, and for a given integer N, define the time step δ t = T/N and the time instants t n = nδ t , n = 0, . . . , N. The value u * k (·, t n ) is denoted by u n+1 k , and let us approximate the time derivative by using the backward Euler scheme, to obtain the following (hereafter, omitting the ' * ' from the local velocity): Next, letû k be a first approximation to u n+1 k , adding the zero −û k +û k to the first term on the left hand side allows us to split the above equation as follows: Furthermore, multiplying Equation (12b) by leads to: By applying gradients on both sides and because u n+1 satisfies the continuity Equation (4), the pressure equation is obtained: Next, by re-using Equation (12b), we formulate the following semi-discrete explicit problem: Given the initial state u 0 k (·), for each n = 0, 1, . . . , N − 1, find the tentative velocityû k , the pressure P n , and update the velocity u n+1 k by solving the following system of semidiscretized equations: The above system is uncoupled, so given u n k , k = 1, 2 we obtainû k from Equation (15a), following the pressure P n from (15b), and finally update u n+1 k with Equation (15c). Now we address its variational formulation, for this, multiply Equations (15a) and (15b) by test functions v, q ∈ V, integrate them by parts, and assuming thatû k , u n k , P n ∈ V, we can formulate Note that the expansion µ ( ∇u n k + u n k ∇ ) is due the variable porosity = (x), and the inner products defined on the boundary Γ result from the integration by parts.
Regarding the Dirichlet boundary conditions (6a,b), they were imposed in a numerical level, so the natural condition n · ∇ u n+1 k = 0 has to be applied on Γ in ∪ Γ w . Meanwhile, to fulfill the condition on the outlet (6c), we split it into a velocity part and a pressure part, n · ∇ u n+1 k = 0 and P n = 0, respectively, enforcing the latter with a generalized condition n · ∇P n = σP n , here σ must be sufficiently large. Finally, the condition for the pressure on walls and the inlet boundary (6d) is applied.
Next, we formulate the finite element spaces associated with triangulation τ h and we define the full discretized problem: Given the initial state u 0 k , k = 1, 2, find the vector (u n+1 1 , u n+1 2 , P n ) ∈ W h × W h × M h that satisfies the follow explicit system for each index n = 0, 1, 2, . . . , N and for all pairs By introducing the nodal basis of spaces W h , M h , it is well known that problem (17) is then equivalent to solving a system of lineal equations built with the contribution matrices of each triangle element in τ h . (17) is the local velocity on the REV u * k at instant t n , and not the average velocity u n k required in the model (7). Therefore, the average velocity has to be updated as follows:

Remark 4. Is important to consider that the unknown in the problem
for each index n = 1, 2, . . . , N.

A Finite Element Approach and Implicit Time Schemes to Solve the Heat-Transfer Model
The model (7) has the classical form of a transport problem, so unlike the Chorin method, we follow the standard treatment of the finite element method and start by performing the spatial discretization. For this, Equation (7a,c) is multiplied by a test function v ∈ V, integraed by parts, and the boundary condition (10a,b) is applied. This leads to the variational formulation of the problem: It must be considered that the above equations, along with Equation (7b), form a nonlinear system, so computing the gradients of the not-yet-defined residuals is essential for solving it. Consider the finite element space associated with triangulation τ h : Therefore, substituting V by V h , taking v ∈ V h , and replacing θ(x, t) by θ(t)w(x) with w ∈ V h leads to: where f a (θ a , θ 0 ) denotes the forcing by convection and radiation (right hand side of (19a)). At this point, several time-marching schemes can be applied to resolve both equations. In this work, we chose a hybrid scheme with a BDF method of order 2 [26] to obtain stability at the first equation and the backward Euler for the last equation. Both schemes are implicit and then unconditionally stable. However, the nature of each equation is non-linear, nonlinear algebraic, and linear, respectively. This motivated us to formulate an implicit method characterized by a time interval between the state variable of each equation and the other unknowns of the system, in order to use the available information.
Considering all of this, it is possible to formulate the complete discretized problem: Given the initial conditions θ 0 a and θ 0 s , find the solutions θ n a , θ n 0 , and θ n s for time instants n = 0, 1, 2, . . . , N, such that it satisfies the following system: The first equation requires two initial values: The first one is the initial condition, and the second one was estimated with a suitable Runge-Kutta method. In this work, we used an implicit Euler scheme. At this stage, it is possible to compute the gradients needed to solve the non-linear Equation (22a,b). Therefore, let R a and R 0 be the residuals defined from Equations (22a,b), and let their gradients (while omitting index n) be expressed as follows: Similarly, by introducing the nodal basis of space V h , a non-linear system can be built with the contribution matrices defined by each term of Equation (22a,c) and vectors from the non-linear vector Equation (22b). The solution of this system produces the values of the state variables in each mesh node and each time instant. Finally, MATLAB ™ software was used to implement the finite element method and, particularly, its command fsolve was used to address the system of non-linear equations. This solver allowed us to supply the gradient (23a,b) and optionally evaluate its accuracy by comparing them with the solver's gradients (estimated with central finite differences). This comparison proved to be effective, with a maximum error of order O(10 −5 ) and O(10 −6 ) for (23a) and (23b), respectively.

Parameters Values and the Urban-Rural Domain Based on the Metropolitan Region of Guadalajara
The domain of study considered in this work is the metropolitan region of Guadalajara in Mexico (shown in Figure 2a). Here, we observed an irregular urban boundary, along with obstacles and hills. The domain Ω was bounded by a rectangle Ω h , surrounding the hills and other small rectangles, and the constructed geometrical domain was meshed by the triangulation τ h (Figure 2b). It contained 5,360 vertices and 10,338 triangular elements such that they were finer on boundary hills and urban regions in order to minimize the consequences of the no-slip boundary condition (6b) and the large temperature gradients expected in a UHI context. Moreover, we assumed that the urban-rural domain had a concentric distribution of buildings, so the center of the city presented a higher concentration of large buildings, in contrast with the suburbs, where smaller and more scattered structures were more common, and with the rural regions, where building concentration was practically nil. Therefore, the notable discrepancy on the model parameter values were intended to agree with our idealized urban-rural domain.
A key point of this work was to consider a variable porosity to characterize the urban and rural regions, instead of the use of urban-limit layouts (and without the need of boundary conditions). Therefore, given the typical porosity value [18,20] for the urban area urb = 0.38 that is associated with high building concentrations, and for the rural area, rural = 0.98 is associated with the lack of building structures. We used a simple Gaussian distribution with variances α x , α y , and centered at the city center (x c , y c ) (see Figure 3b) to distribute them across the domain. Furthermore, in order to define the urban limits, we modeled the city as a circle with radius r, which we used to truncate the Gaussian distribution at the points where it reaches the city limits (see black and magenta lines on Figures 2b and 3a). The porosity was not the only parameter related to the concentration of building blocks, as there were other factors, such as soil density (ρ urb = 2.11 × 10 3 , ρ rural = 8.4 × 10 2 ), soil roughness (z 0,urb = 7, z 0,rural = 1), Bowen radius (β urb = 5, β rural = 0.5), and others that were defined with a Gaussian distribution. Nevertheless, there were dimensionless parameters, such as emissivity (e urb = 0.84, e rural = 0.96), and albedo (a urb = 0.27, a rural = 0.16), which were not dependent on the concentration of buildings. Therefore, they were defined with a radial distribution (see Figure 3b). Lastly, we used conservative values for air and soil layer thickness d a = 2 m and d s = 1 m. For a complete list of the model parameters, and the parameters associated with either urban or rural regions, see Table A1.

Remark 5.
For the sake of simplicity, we used an idealized urban-limit layout. More realistic urban layouts could be built with the curve fitting of a finite number of points on the real city limits. Nevertheless, the described methodology applies in such realistic cases.

A 24 h Simulation of the UHI under Ideal Conditions
In our first test, we considered a simulation of T = 24 h under ideal conditions: no wind, a typical day-long radiation (Figure 2c), and no clouds, which was carried out with the aim of reproducing the temperature inversion described in [3,11]: a warm city during the day that transitioned into a cold city at night. Results are shown in Figure 4, where the spatial distribution of temperature θ a at 8, 10, 12, 14, 16, 18, and 20 h of the day is displayed. We observed not only the expected transition from a warmer urban process (consistent with radiation increase) to a colder process (consistent with radiation decrease). Both were in contrast with their rural surroundings, but they also had different heat levels, with differences of 2-3 degrees in the city (Figure 4c). This was the intended inverse temperature effect. Remark 6. It could be accurate to say that [11] observed the temperature inversion on parks surrounded by an urban area, that is, parks turning from a cold island during the day to a warmer island at night, with respect to their urban surroundings. Given that parks are conformed by plants, trees, and green soil, that have similar thermal properties as the rural surroundings of the city, this allowed us to assume that this phenomenon holds for a city surrounded by a rural (green) area. Therefore, the city transitioned from a heat island during the day to a cold island at night, with respect its rural surroundings.

Influence of the Wind and the Need for Numerical Stabilization
With respect to the wind field, we solved the Darcy-Forchheimer-Brinkman model by forcing it with the boundary condition values g 1 = 0.25 m s −1 , and g 2 = −0.25 m s −1 , and then with an inlet wind gust in a southwest direction. Therefore, with enough time, the model achieved a stationary state which was used as the reference wind (see Figure 5), having maximum magnitude of 1.2 m s −1 . To augment the wind intensity, we used a scalar parameter η that multiplied its intensity by 1, 2, 3, 4 times. The model is highly convective due to the low density of the air and the intensity of the wind field, and this combination tends to destroy the stability of the numerical solution through spurious oscillations. This is shown in Figure 6a, where the air temperature oscillations were present along the wind direction when η = 4 and the wind blew for 3 h. To address this problem, several techniques could have been used. We chose the streamline diffusion stabilizer [28], and as shown in Figure 6b, the oscillations were minimal. The streamline diffusion stabilizer added diffusion, but only in the direction of the wind. By catching the temperature changes using an inner product between wind and θ gradient, it was implemented in the variational formulation level by adding a suitable term to Equation (19a): , v + u · ∇θ a , v + µ a ∇θ a , ∇v + ζ u · ∇θ a , u · ∇v = f a (θ a , θ 0 ), v , where ζ > 0 is a parameter to modulate the stabilization. Certainly, the stabilizer is artificial, and the new variational formulation is not consistent with the strong formulation of the model, but these types of adjustments are not exclusive to the Galerkin method. For example, discontinuous Galerkin and finite volumes in heat transfer and fluid dynamics both use WENO stabilizers, which are implemented to limit the gradients between the elements of the mesh [32][33][34]. Therefore, in view of the results shown in Figure 6, future numerical experiments should use the streamline diffusion stabilizer.
Next, we analyzed the influence of different wind intensities η = 1, 2, 3, 4 blowing from hour 10 to hour 13 in a day. The results are shown in Figure 7, where, as expected, the wind carried the hot mass of air to the southeast and, thereby, benefiting the northwest city inhabitants (mainly those living close to the city limits) with fresh air but disadvantaging those living in the southeast and outside of the city. The differences in the air temperatures could reach 2 • K. We had to consider that radiation between the hours 10 and 13 increased, thereby continuing to heat the system. The city center did not show significant changes and maintained higher temperatures. Below, we show the changes during the afternoon. Our last test demonstrated an interesting feature of the UHI, as it occurred when the radiation was low enough to stop the warming of the system, but at the same time, it was still subject to intense wind. The model output showed the residual mass of warmer air from the city continued to move through the countryside, exhausting its energy by interchanging it with the colder surface, until it reached the environmental temperature. At specific instants, warmer air masses had a temperature as high as 3 • K higher than the rural surroundings, which was unrelated to this dynamic. For an example, see the southeast zone (Figure 8b,c). In contrast, the UHI was broken in the northwest part of the city, where the temperature contours overlapping with the city circumference boundary were no longer present. We noted that the similarities between the city circumference with temperature contours, as shown in Figure 8a, were no longer showing in the northwest section, as shown in Figure 8b, but it returned in Figure 8c with colder temperatures.

Discussion
In this work, a numerical study of the well-known phenomenon, urban heat island (UHI), was developed using key contributions and could be compared with other works. The characterization of the rural and urban regions using a spatially varying porosity defined by a Gaussian distribution was a different approach, with respect to the urbanmonolithic type, with constant parameters or buildings distribution depending on the porosity used in the previous works ( [12,[16][17][18]20] and the references therein). It was interesting that the variability of the parameters enabled us to obtain a notable contrast, not only between the city and its rural surroundings, but also in the city itself, with temperature discrepancies of a few degrees. However, the porosity we defined had consequences on the model, its variational formulation, and the finite element approach, which have been explained in detail throughout the paper.
With respect to the CFD and heat-transfer models and their numerical solutions, a simple adjustment in the Darcy-Forchheimer-Brinkman equations after rewriting them with respect to the local velocity u * k , instead of the average fluid velocity u k , ensured they became classic equations of fluid dynamics with porosity-dependent forcing. Furthermore, the model for heat transfer was written based on equations used in remote sensing [23] but was adapted to the energy interchange in porous media [21]. Both allowed us to reach the proposed goals: temperature inversion and warm air transport.
From our perspective, that both models were resolved with classical numerical schemes with great accuracy was remarkable and should enable other researchers to reproduce and amplify our results. This is in contrast with more complicated methodologies, such as discontinuous finite element and finite volume methods and their limiters.
With respect to expanding our work, we inferred that our model would be able to simulate the effects of the strategies in order to mitigate the intensity of the UHI, which was evidenced in recently published reviews [15,35]: great and median urban parks [35,36], the use of white-green roofs [37], trees with dense vegetation [15,38], and the use of concrete instead of asphalt [39]. However, those strategies must be carried out in a large enough area to influence the numerical results.
Last, but not least, this paper presented challenges to be considered in future research, such as a sensitivity analysis of the influence of the model parameters on the numerical outputs; a vertical extension of the model that considers the turbulent transport; and improvements to the layout of the urban limits. Regarding the latter, there are at least two possible alternatives, using boundary conditions in an urban layout (more complicated mathematical formulation and computational effort) or defining the urban limit layout by truncating a continuous function, as we did in this work.