Multiscale Finite Element Method for Urban Canopies in Climate Models

By Heena Patel

A climate model is a mathematical representation of the climate system that is based on physical, biological, and chemical laws. The equations that are derived from these laws are very complex and thus require numerical solutions. Climate models therefore provide a discrete solution in space and time — namely, averages for specific times over regions whose sizes depend on the model’s resolution. Some models may only provide globally or zonally averaged values, while others may have numerical grids with spatial resolution that are smaller than 100 kilometers. The time step might range from a few minutes to several years. Even with the highest resolution, the numerical grid is too coarse to represent small-scale processes like turbulence in the atmospheric and oceanic boundary layers, the circulation’s interaction with small-scale topography features, thunderstorms, canopy microphysics processes, and so forth. Furthermore, many processes are still not understood sufficiently enough to enable detailed modeling of their behavior.

In order to account for the influence of small-scale processes on larger scales, researchers must design parameterizations that are based on empirical evidence and/or theoretical arguments for processes that are not explicitly included. These parameterizations only reproduce first-order effects and are not valid for all possible conditions; as such, they are often a significant source of uncertainty in models [4]. Consequently, the development of a mathematically consistent method for the transfer of information from the subgrid scale to the coarse grid is necessary. Multiscale numerical modeling offers a promising mathematical framework for achieving this goal.

The multiscale finite element method (MsFEM) consists of two major ingredients: (i) multiscale basis functions that capture the solution’s fine-scale features and (ii) a global numerical formulation that couples these multiscale basis functions. We must incorporate important multiscale features of the solution into the localized basis functions, which contain information about the scales that are smaller and larger than the local numerical scale that is defined by the basis functions. In particular, we need to incorporate the features of the solution that can be localized and use additional basis functions to capture information about the features that are included separately in the coarse space. A global formulation couples these basis functions to provide an accurate approximation of the solution [3].

This type of research aims to study and implement MsFEM in urban climate simulations and consider the advection-diffusion equation to observe various changes in flow parameters—such as velocity, temperature, and pollutants—thereby extending the work in [7, 8, 9]. Here we use deal.II, a C++ program library that targets the computational solution of partial differential equations [1]. This advanced mathematical library incorporates massively parallel high-performance computing features for various applications and has several advantages, including user control of mathematical implementation, availability of different preconditioner solvers, and detailed documentation. Using deal.II, we implemented a parallelization algorithm for semi-Lagrangian meshes on different processors to support the application of MsFEM to transport-dominated problems.

Scientists develop various parameterization schemes in order to implement the canopy effect in numerical weather prediction or climate models. One study explores the numerical aspects of physical parameterization within the context of the European Centre for Medium-Range Weather Forecasts Integrated Forecasting System [2]. In this case, researchers must implement clouds and land cover into the equation by adding source terms. Additionally, various processes have different numerical problems. In the context of numerical design, one should thoroughly understand the physics of the problem and the interaction between different processes. Another study found that more accurate canopy representations incur a higher computational cost when different urban canopy models are coupled with weather forecast models [5]. Understanding the effect of canopy at various scales therefore becomes necessary. Researchers have previously implemented canopies in the advection-diffusion equation via a function of diffusion coefficients [6, 10].

Here, we represent buildings as obstructions to demonstrate the difference in diffusivity and the effect of tracer transport (see Figure 2). The setup mimics a real situation with a single building under wind and concentration conditions. We aim to implement a simple stationary obstacle in the advection-diffusion equation and verify the results with a wind tunnel experiment at the University of Hamburg’s Environmental Wind Tunnel Laboratory. These efforts represent a simple step towards MsFEM application, from the microscale (where the canopy lies) towards larger scales.

Consider the following advection-diffusion equation with periodic boundaries for two and three dimensions:

\[\partial _t u + c_\delta \cdot \nabla u = \nabla \cdot (A_\epsilon \nabla u) + f\ \ \ \textrm{in}\ \ \ \textrm{T}^\textrm{d} \times [0,\textrm{T}]\sim [0,1]^\textrm{d} \times [0,\textrm{T}]\tag1\]

\[u(x,0) = u_0 (x)\]

and

\[\partial _t u + \nabla \cdot (c_\delta u) = \nabla \cdot (A_\epsilon \nabla u) + f\ \ \ \textrm{in}\ \ \ \textrm{T}^\textrm{d} \times [0,\textrm{T}] \sim [0,1]^\textrm{d} \times [0,\textrm{T}]\tag2\]

\[u(x,0) = u_0 (x),\]

where \(c_\delta\) is the background velocity. \(A_\epsilon (x,t)\) is a matrix-valued diffusivity, and \(f\) and \(u_0\) are smooth external forcing and initial conditions. The equations can be transformed into a discrete set of statements of the form \(u^{n+1} = F(u^n)\), where \(F\) is some functional that discretizes equations \((1)\) or \((2)\) and \(u^n\) is the discrete representation of the unknown solution \(u\) at time steps \(t^n\). The algorithm for the finite element method (FEM) is thus as follows:

1. Set up mesh
2. Set up system and constraints (set \(n=0\))
3. While \(t^{n+1} \le T\), assemble the system, solve for \(u^{n+1}\), and set \(n=n+1\).

Accordingly, the algorithm for MsFEM takes the following steps:

1. Set up coarse mesh
2. Set up system and constraints
3. Compute Ms-basis at \(t=0\) (set \(n=0\))
4. While \(t^{n+1} \le T\), compute basis at \(t^{n+1}\), assemble system, solve for \(u^{n+1}\), and set \(n=n+1\).

Figure 2 depicts the reference solution for a single building with a height of 30 meters. In Figure 3, we replicate this setup in a wind tunnel experiment; a black block represents the building and red Lego blocks add an element of roughness.

We parameterize the building with a diffusion coefficient maximum inside the structure and an outside linear function. The velocity is set to zero inside the building, grows exponentially in a small transition layer near the building to a preset profile in the horizontal direction, and grows logarithmically above the building in the vertical direction, as is expected for a boundary layer profile. Figure 4 illustrates low-resolution MsFEM’s ability to capture subgrid processes more accurately than standard FEM, though less accurately than high-resolution FEM which requires many more computations.

In conclusion, the canopy exists at a subgrid scale in climate models but plays an influential role in the modulation of local—but large-scale—climate. Canopy modeling has become increasingly significant as global temperatures continue to rise above preexisting limits. While numerous studies have taken place in this area, they incur high computational costs. In order to account for canopy effects, models include canopy parametrization as a source term in flow or transport equations. Researchers typically add a source or sink term to the equation in order to represent canopies — i.e., buildings are considered to be solid and are surrounded by the fluid domain in the mesh, based on the different boundary conditions that are applied for the fluid and solid domains. They then solve the system using either terrain-following coordinates or the immersed boundary method. Mesoscale models also account for external forcing of data from climate simulation. All of these methods are computationally expensive when the models are run at a high resolution.

Our approach simplifies this complexity by representing the obstacles in the canopy with a (large) diffusion coefficient and employing a MsFEM discretization for the subgrid-scale structure of the canopy layer. We tested our approach by conducting wind tunnel experiments with simple one-building setups and comparing the data with numerical solutions. While our experiments are idealized cases that demonstrate the general feasibility of this new method, a generalization promises to enable new ways to solve the so-called upscaling problem, in which small-scale features—which are not representable numerically—must be represented in large-scale dynamics. Thus, we established the potential for future implementation in global climate models.

---

Heena Patel presented this research during a contributed presentation at the 2022 SIAM Annual Meeting (AN22), which took place this July in Pittsburgh, Pa. She received funding to attend AN22 through a SIAM Student Travel Award. To learn more about Student Travel Awards and submit an application, visit the online page. 

Acknowledgments: This work was conducted under the supervision of Konrad Simon and Jörn Behrens at the University of Hamburg. The author extends a special thanks to Sylvio Freitas, who conducted the wind tunnel experiment, as well as Frank Harms and Bernd Leitl from the University of Hamburg’s Environmental Wind Tunnel Laboratory. This research is funded by Cluster of Excellence “Climate, Climatic Change, and Society” (CLICCS) Project A3 - Canopies in the Earth System. 

References
[1] Alzetta, G., Arndt, D., Bangerth, W., Boddu, V., Brands, B., Davydov, D., … Wells, D. (2021). The deal.II library, version 9.3. J. Numer. Math., 26(4), 173-184.
[2] Beljaars, A., Balsamo, G., Bechtold, P., Bozzo, A., Forbes, R., Hogan, R.J., … Wedi, N. (2018). The numerics of physical parametrization in the ECMWF model. Front. Earth Sci., 6.
[3] Efendiev, Y., & Galvis, J. (2012). Coarse-grid multiscale model reduction techniques for flows in heterogeneous media and applications. In I. Graham, T. Hou, O. Lakkis, & R. Scheichl (Eds.), Numerical analysis of multiscale problems (pp. 97-125). Lecture notes in computational science and engineering (Vol. 83). Heidelberg: Springer
[4] Goosse, H., Barriat, P.Y., Lefebvre, W., Loutre, M.F., & Zunz, V. (2008-2010). Introduction to climate dynamics and climate modeling. Louvain-la-Neuve, Belgium: Université Catholique de Louvain. Retrieved from http://www.climate.be/textbook.
[5] Jandaghian, Z., & Berardi, U. (2019). Proper choice of urban canopy model for climate simulations. In 16th IBPSA international building performance simulation conference and exhibition. Rome, Italy: International Building Performance Simulation Association.
[6] Kanda, I., Yamao, Y., Ohara, T., & Uehara, K. (2013). An urban atmospheric diffusion model for traffic-related emission based on mass-conservation and advection-diffusion equations. Environ. Model. Assess., 18, 221-248.
[7] Simon, K., & Behrens, J. (2019). A semi-Lagrangian multiscale framework for advection-dominant problems. In Computational Science – ICCS 2019 (pp. 393-409). Faro, Portugal: Springer.
[8] Simon, K., & Behrens, J. (2020). Multiscale finite elements for transient advection-diffusion equations through advection-induced coordinates. Multiscale Model. Simul., 18(2), 543-571.
[9] Simon, K., & Behrens, J. (2021). Semi-Lagrangian subgrid reconstruction for advection-dominant multiscale problems with rough data. J. Sci. Comput., 87, 49.
[10] Park, Y.-S., & Paw U, K.T. (2004). Numerical estimations of horizontal advection inside canopies. J. Appl. Meteorol., 43(10), 1530-1538.

Heena Patel is a Ph.D. student in the Numerical Methods in Geosciences group at the University of Ham-burg. Her focus is the implementation of the multiscale finite element method for application to urban canopies through high-performance computing.