QES-Winds

Introduction

A new dispersion modeling system based on the well-used FORTRAN-based QUIC (Quick Urban and Industrial Complex) dispersion modeling system originally developed by the University of Utah and Los Alamos National Laboratory [1], has been under development as collaboration between the University of Utah and the University of Minnesota, Duluth. Quick Environmental Simulation (QES) is a microclimate simulation platform for computing 3D environmental scalars in urban areas and over complex terrain.

The fast-response three-dimensional diagnostic wind model written in C++, QES-Winds, is a rapid mass conserving wind-field solver. QES-Winds utilizes the concept of dynamic parallelism in CUDA to substantially accelerate wind simulations. Figure below shows a high-level flowchart for QES-Winds code.

Flowchart for the QES-Winds wind solver

Mass Consistent Solver

Staggered Grid

QES-Winds discretizes the computational domain using a staggered grid as shown in figure below.

Staggered grid representation of the domain and location of each variable.

The velocity components (\(u\), \(v\) and \(w\)) are face-centered values and Lagrange multipliers (\(\lambda\)), divergence of the initial wind field (\(R\)) and solver coefficients (\(e\), \(f\), \(g\), \(h\), \(m\) and \(n\)) are cell-centered variables. Because of nature of the finite difference method (depending on neighboring cells values), the first and last cells in \(x\) and \(y\) directions and the last cell in \(z\) direction, are not updated in the solving process (their velocity remains as the same as the initial velocity field). For the same reason, there is a layer of ghost cells under the bottom surface to make velocity calculation in the first layer above the bottom surface possible. The values of the Lagrange multipliers for the ghost cells are set to the ones for the layer above the bottom surface to create a zero gradient for the Lagrange multipliers (boundary condition) as well as providing the neighboring cell for the finite difference method.

Poisson Equation

QES-Winds have mass conserving wind field solvers that rapidly compute wind fields using a variational method rather than slower yet more physics based solvers that include conservation of momentum [3]. While the QES-Winds method uses reduced order physics in the numerical solution of urban flow problems, the solutions are rapid and compare quite well higher order physics-based models in both idealized [4] and realistic urban cities [5]. The method minimizes the difference between an initial wind field that is specified using empirical parameterizations [6] and the final wind fields. The empirical parameterizations account for complex wind fields around buildings such as wake cavities downstream of a building. To obtain a quasi-time-averaged velocity field, QES-Winds uses a variational analysis technique [6]. This method requires the solution of a Poisson equation for Lagrange multipliers, \(\lambda\) in the following form:

\[\label{poisson} \frac{\partial^2\lambda}{\partial x^2} + \frac{\partial^2\lambda}{\partial y^2} + (\frac{\alpha_1}{\alpha_2})^2\: \frac{\partial^2\lambda}{\partial z^2} = R\]

Where R is divergence of the initial wind field and is defined as:

\[\label{divergence} R = -2\,\alpha_1^2\,\Bigg[\frac{u_{i+1/2}^0-u_{i-1/2}^0}{\Delta x} + \frac{v_{j+1/2}^0-v_{j-1/2}^0}{\Delta y} + \frac{w_{k+1/2}^0-w_{k-1/2}^0}{\Delta z}\Bigg]\]

The final velocity field is updated using Euler-Lagrange equations:

\[\label{eu-lag1} u = u^0 + \frac{1}{2\,\alpha_1^2\,\Delta x}\,[\lambda_{i+1\,,j,\,k}-\lambda_{i,\,j,\,k}]\]

\[\label{eu-lag2} v = v^0 + \frac{1}{2\,\alpha_1^2\,\Delta y}\,[\lambda_{i,\,j+1,\,k}-\lambda_{i,\,j,\,k}]\]

\[\label{eu-lag3} w = w^0 + \frac{1}{2\,\alpha_2^2\,\Delta z}\,[\lambda_{i,\,j,\,k+1}-\lambda_{i,\,j,\,k}]\]

The Poisson equation is solved using the Successive Over-Relaxation (SOR) method which is a variant of Gauss-Seidel method with faster convergence. Applying the SOR scheme to the Poisson equation for Lagrange multipliers results in:

\[\begin{split}\label{SOR} \begin{split} \lambda_{i,\,j,\,k} & = \frac{\omega\Bigg[(\Delta x)^2 R_{i,\,j,\,k}+e\,\lambda_{i+1}+f\, \lambda_{i-1}+A(g\,\lambda_{j+1}+h\, \lambda_{j-1}) + B(m\,\lambda_{k+1}+n\, \lambda_{k-1})\Bigg]}{e+f+g+h+m+n}\\ & +(1-\omega)\lambda_{i,\,j,\,k} \end{split}\end{split}\]

Where e,f,g,h,m,n are boundary condition coefficients and A and B are domain constants. \(\omega = 1.78\) is the SOR relaxation factor. The boundary condition for solid surfaces is (\(\frac{\partial \lambda}{\partial n}=0\)) and for inlet/outlet surfaces it is \(\lambda=0\).

Solver Types

QES-Winds has four options for solving the SOR equation discussed above, the first option is to solve the equation on the CPU and the rest use the GPU for computations. The GPU solvers are called: the dynamic parallel, the global memory and the shared memory. The CPU solver is quite rapid, but slow in comparison to the GPU solvers since it is a serial solver and does not have parallel computing capabilities, especially for large domains. For more information regarding different types of solvers available in QES-Winds, read [2].

Parameter Files

The XML parameter file has the following structure, with the XML elements corresponding different section of the model. Each of them are presented in the sections below (expect <turbParams> which is presented in QES-Turb).

<QESWindsParameters>
    <simulationParameters>
        <!-- HERE COMES THE SIMULATION PARAMETERS -->
    </simulationParameters>

    <metParams>
        <!-- HERE COMES THE MET PARAMETERS -->
    </metParams>

    <buildingsParams>
        <!-- HERE COMES THE BUILDING PARAMETERS -->
    </buildingsParams>

    <vegetationParams>
        <!-- HERE COMES THE VEGETATION PARAMETERS -->
    </vegetationParams>

    <turbParams>
        <!-- HERE COMES THE TURBULENCE PARAMETERS -->
    </turbParams>

    <fileOptions>
        <!-- HERE COMES THE FILE PARAMETERS -->
    </fileOptions>
</QESWindsParameters>

QES-Winds Domain (simulationParameters)

The first step in every computational code or package is to define the computational domain. The user can define the domain by specifying the number of cells in \(x\), \(y\) and \(z\) directions as well as the cell size in each direction in the input file (XML file).

Basic Parameters

The domain information (number of cells and cell size) are defined under the <simulationParameters> part of the XML file. Following is an example of a domain with \(2\) km by \(2\) km by \(200\) m and resolution of \(2\) m by \(2\) m by \(2\) m:

<simulationParameters>
    <!-- Number of cells in x,y and z directions-->
    <domain> 1000 1000 100 </domain>
    <!-- Mesh resolution (meters)-->
    <cellSize> 2.0 2.0 2.0 </cellSize>
    <!-- vertical stretching (0-uniform grid (default), 1-costum grid)-->
    <verticalStretching> 0 </verticalStretching>
    <!-- Number of time steps-->
    <totalTimeIncrements> 1 </totalTimeIncrements>

    <!-- ... -->
</simulationParameters>

Note, <verticalStretching> is currently under development. The last parameter lets the user define the number of time instances that need to run, assuming the correct number of sensor time steps are define in <metParams>. This parameter can be set to 0 to let the program decide how many time instances need to be run based on input parameters.

Halo Region

If a solid element (building or terrain) overlaps with the QES domain boundaries, QES-Winds cannot model the wind field around the element correctly. In order to prevent this phenomenon, the user can add buffer zones to the sides of the domain when a terrain file or an ESRI shapefile is read into the code. Figure below represents how the halo region is added to the domain around a Digital Elevation Model (DEM) or a shapefile.

Representation of halo region around the domain.

In order to define length of the halo zone in \(x\) and \(y\) direction, the user can use <halox> and <haloy> under <simulationParameters>. When the halo zone is defined, the length of the domain (\(nx*dx\)) and (\(ny*dy\)), should be greater than or equal to length of the DEM or shapefile in each direction plus twice the length of the halo in \(x\) and \(y\) directions, respectively.

<simulationParameters>
    <!-- Halo region added to x-direction of domain (at the beginning and the end of domain) (meters)-->
    <halo_x> 20.0 </halo_x>
    <!-- Halo region added to y-direction of domain (at the beginning and the end of domain) (meters)-->
    <halo_y> 30.0 </halo_y>

    <!-- ... -->
</simulationParameters>

Digital Elevation Model (DEM)

The current version of QES-Winds has been written to allow commonly available terrain and building geometry datasets to be used for simulations. In this section, various input file formats for QES-Winds will be covered.

Terrain Features

Using the Geospatial Data Abstraction Library (GDAL; https://www.gdal.org), we are able to load geo-referenced datasets of terrain so that the simulations can include the effects of hills, valleys, and mountains. In the current version of the code, we can load Digital Elevation Model (DEM) files for different physical locations.

Using the Digital Elevation Model (DEM) file loaders in our code base, we have loaded and tested multiple different terrain data sets. As a first test, we loaded a DEM of Askervein Hill. This is an isolated hill in Scotland where field experiments have been conducted and data for testing and evaluation exists ([7, 8]). The simulation with Askervein Hill was run without any complex terrain flow parameterizations. The Askervein Hill dataset is \(6023.43\) m by \(6023.43\) m. The hill height is approximately \(124\) m tall. Figure below indicates the cell type contour for the Askervin hill test case in a vertical plane at \(y = 3000\) m (part (a)), and a horizontal plane at \(z=20\) m (part (b)). These plots show the ability of QES-Winds to read in and process DEM files. The cell type value \(1\) (blue) represents the air cells while value \(2\) (red) indicates the terrain cells.

Cell type contour for the Askervin hill test case in a (a) vertical plane at \(y=3000\) m, (b) horizontal plane at \(z=20\) m. The cell type value \(1\) (blue) represents the air cells while value \(2\) (red) indicates the terrain cells.

The user can define the address to the DEM using <DEM> variable under the <simulationParameters> part in the XML file:

<simulationParameters>
    <!-- Address to DEM location-->
    <DEM>../scratch/DEM/askervein.tif</DEM>

    <!-- ... -->
</simulationParameters>

Process Part of DEM

In some cases, user wants to load a giant DEM but only process part of the file. This is possible in QES-Winds by defining the origin of QES domain inside the DEM borders and the size of the QES domain. Figure below shows a schematic of how the QES domain can be defined inside a DEM file and only process that part.

Schematic of how the QES domain can be defined inside a DEM file and only process that part.

There are two options to determine the location of the origin of QES domain inside the DEM borders:

1. Specifying the distance of the QES origin with respect to bottom left corner of the DEM file. This can be done by setting the value of <originFlag> to \(0\) and defining distances (in meters) in \(x\) and \(y\) directions using <DEMDistancex> and <DEMDistancey>, respectively.

   <simulationParameters>
       <!-- Origin flag (0- DEM coordinates (default), 1- UTM coordinates) -->
       <originFlag> 0 </originFlag>
       <!-- x component (m) of origin in DEM coordinates (if originFlag = 0) -->
       <DEMDistancex> 1000.0 </DEMDistancex>
       <!-- y component (m) of origin in DEM coordinates (if originFlag = 0) -->
       <DEMDistancey> 1000.0 </DEMDistancey>
   
       <!-- ... -->
   </simulationParameters>
   

2. Defining the location of the QES domain origin in the Universal Transverse Mercator (UTM) coordinates by setting the value of <originFlag> to \(1\) and determining <UTMx> and <UTMy> of the origin in \(x\) and \(y\) directions, respectively.

   <simulationParameters>
       <!-- Origin flag (0- DEM coordinates (default), 1- UTM coordinates) -->
       <originFlag> 1 </originFlag>
       <!-- x component (m) of origin in UTM DEM coordinates (if originFlag = 1)-->
       <UTMx> 595469.6122881 </UTMx>
       <!-- y component (m) of origin in UTM DEM coordinates (if originFlag = 1)-->
       <UTMy> 6336281.9538635 </UTMy>
   
       <!-- ... -->
   </simulationParameters>
   

Initial Wind Field (metParams)

QES-Winds can read a single or multiple sensors for a specific test case. In this context, sensor means the velocity magnitude and direction at a single point or a single velocity profile to initialize the wind field. If there is only the wind velocity and direction at a single point, the user should specify what type of velocity profile they want to build from the measurement. There are three options available for the type of profile:

1. a logarithmic profile [9]:

   \[\label{eq:log_law} u_{log}(z) = u_{ref}\cdot\frac{\ln(z/z_0)}{\ln(z_{ref}/z_0)}\]

   where \(u_{ref}\) is the measured velocity at measured height \(z_{ref}\), \(z_0\) is the surface roughness.

2. a power law profile [9]:

   \[\label{eq:power_law} u_{pow}(z) = u_{ref}\cdot(z/z_{ref})^{p}\]

   where \(u_{ref}\) is the measured velocity at measured height \(z_{ref}\), \(p\) is the power law exponent.

3. an urban canopy profile [9, 10]:

   \[\begin{split}\label{eq:urban_canopy_low} u_{uc}(z)=\begin{cases} u_H\cdot\exp(\alpha(\frac{z}{H}-1)) & \text{if } z\leq H\\ \frac{u_*}{\kappa}\cdot \ln(\frac{z-d}{z_0}) & \text{if } z > H. \end{cases}\end{split}\]

   Here the profile is separated into two parts: below and above the canopy height \(H\). The method uses the reference velocity \(u_{ref}\) is the measured velocity at measured height \(z_{ref}\) and the attenuation coefficient \(\alpha\) to calculate, via bisection, the zero plane displacement height \(d\), the velocity at the canopy top \(u_H\), and the friction velocity \(u_*\). The lower portion of the urban canopy profile is calculated as an exponential reaching the velocity at the canopy top \(u_H\), and the upper portion of the urban canopy is a logarithmic profile using the computed \(u_*\) and \(d\), with \(\kappa \sim 0.4\) the von Karman constant.

Figure below shows velocity profiles of used for the initial velocity field using each of the scheme presented above.

Velocity profiles created using (1) logarithmic profile, (2) power law profile, and (3) urban canopy profile.

If there is only one sensor available in the computational domain, the code will extend the profile for that sensor uniformly to the whole domain. On the occasion of multiple sensors, QES-Winds utilizes a two-dimensional Barnes interpolation scheme [11, 12] to interpolate velocity components at each cell height of the domain based on the weighted distance from each sensor.

XML Setup

There are two options available for defining sensor information:

1. The user can define all information required for creating a sensor by using the <sensor> variable inside the <metParams> section of the XML file. An example of this section is presented below:

   <!-- Meteorological parameters -->
   <metParams>
       <!-- Distribution of surface roughness for domain (0-uniform (default), 1-custom -->
       <z0_domain_flag> 0 </z0_domain_flag>
   
       <!-- Define a sensor -->
       <sensor>
           <!-- Sensor site coordinate system (1=QES (default), 2=UTM, 3=Lat/Lon) -->
           <site_coord_flag> 1 </site_coord_flag>
           <site_xcoord> 100.0  </site_xcoord>
           <site_ycoord> 140.0 </site_ycoord>
   
           <!-- Start of timestep information for a sensor -->
           <timeSeries>
               <timeStamp>2003-07-23T23:00:00</timeStamp>
               <boundaryLayerFlag> 1 </boundaryLayerFlag>
               <siteZ0> 0.1 </siteZ0>
               <reciprocal> 0.0 </reciprocal>
               <height> 10.0 </height>
               <speed> 5.0 </speed>
               <direction> 270.0 </direction>
           </timeSeries>
       </sensor>
   </metParams>
   

   In that case, multiple sensors can be defined in this section, each added using <sensor>…</sensor>.

2. The user can put all the sensor information in a separate XML file and define the address to the location of the sensor file using the <sensorName> variable.

   <metParams>
       <!-- Distribution of surface roughness for domain (0-uniform (default), 1-custom -->
       <z0_domain_flag> 0 </z0_domain_flag>
       <!-- Name of the sensor file with information for the sensor included -->
       <sensorName>../data/InputFiles/sensor.xml</sensorName>
   </metParams>
   

   An example of the sensor file is presented below:

   <sensor>
       <!-- Sensor site coordinate system (1=QES (default), 2=UTM, 3=Lat/Lon) -->
       <site_coord_flag> 1 </site_coord_flag>
       <site_xcoord> 590.0  </site_xcoord>
       <site_ycoord> 1.0 </site_ycoord>
   
       <!-- Start of timestep information for a sensor -->
       <timeSeries>
           <timeStamp>2003-07-23T23:00:00</timeStamp>
           <boundaryLayerFlag> 3 </boundaryLayerFlag>
           <siteZ0> 0.3 </siteZ0>
           <reciprocal> 0.0 </reciprocal>
           <height> 17.23 </height>
           <speed> 5.15 </speed>
           <direction> 153.83 </direction>
           <canopyHeight> 10.1 </canopyHeight>
           <attenuationCoefficient> 2.41 </attenuationCoefficient>
       </timeSeries>
   
       <!-- ... -->
   </sensor>
   

   In that case, multiple sensors can be defined each in their separate XML file.

The first part of the sensor information is the location of the sensor in domain. There are three options for it:

1. define the location in local coordinates of the QES domain.

   <sensor>
       <!-- Sensor site coordinate system (1=QES (default), 2=UTM, 3=Lat/Lon) -->
       <site_coord_flag> 1 </site_coord_flag>
       <!-- x component of site location in QES domain (m) (if site_coord_flag = 1) -->
       <site_xcoord> 1.0  </site_xcoord>
       <!-- y component of site location in QES domain (m) (if site_coord_flag = 1)-->
       <site_ycoord> 1.0 </site_ycoord>
   
       <!-- ... -->
   </sensor>
   

2. The user can define the location in the UTM coordinates.

   <sensor>
       <!-- Sensor site coordinate system (1=QES (default), 2=UTM, 3=Lat/Lon) -->
       <site_coord_flag> 2 </site_coord_flag>
       <!-- x components of site coordinate in UTM -->
       <site_UTM_x> 634175 </site_UTM_x>
       <!-- y components of site coordinate in UTM-->
       <site_UTM_y> 3925362 </site_UTM_y>
       <!-- UTM zone of the sensor site -->
       <site_UTM_zone> 14 </site_UTM_zone>
   
       <!-- ... -->
   </sensor>
   

   In this case, user also needs to define the origin of computational domain in the UTM coordinates.

   <simulationParameters>
       <!-- x component (m) of origin in UTM -->
       <UTMx> 634173 </UTMx>
       <!-- y component (m) of origin in UTM -->
       <UTMy> 3925360 </UTMy>
       <!-- UTM zone that domain located -->
       <UTMZone> 14 </UTMZone>
   </simulationParameters>
   

3. The user can define the location in Latitude and Longitude coordinates.

   <sensor>
       <!-- Sensor site coordinate system (1=QES (default), 2=UTM, 3=Lat/Lon) -->
       <site_coord_flag> 3 </site_coord_flag>
       <!-- x components of site coordinate in Latitude -->
       <site_lat> 35.46270 </site_lat>
       <!-- y components of site coordinate in Longitude -->
       <site_lat> -97.52130 </site_lat>
   
       <!-- ... -->
   </sensor>
   

   In this case, user also needs to define the origin of computational domain in the UTM coordinates.

   <simulationParameters>
       <!-- x component (m) of origin in UTM -->
       <UTMx> 634173 </UTMx>
       <!-- y component (m) of origin in UTM -->
       <UTMy> 3925360 </UTMy>
       <!-- UTM zone that domain located -->
       <UTMZone> 14 </UTMZone>
   </simulationParameters>
   

The second part of sensor definition is choosing type of profile for different time steps, if applicable. The <timeSeries> variable is designed to define type of sensor profile in the sensor section for several time steps. There are four options for the input profile in QES-Winds:

1. Logarithmic velocity profile:

   <!-- Start of timestep informastion for a sensor -->
   <timeSeries>
       <!-- time of the current data input -->
       <timeStamp>2003-07-23T23:00:00</timeStamp>
       <!-- Site BL flag (1-log (default), 2-exp, 3-urban canopy, 4-data entry) -->
       <boundaryLayerFlag> 1 </boundaryLayerFlag>
       <!-- Site surface roughness -->
       <siteZ0> 0.1 </siteZ0>
       <!-- Reciprocal Monin-Obukhov Length (1/m) -->
       <reciprocal> 0.0 </reciprocal>
       <!-- Height of the sensor -->
       <height> 20.0 </height>
       <!-- Measured speed at the sensor height -->
       <speed> 5.0 </speed>
       <!-- Wind direction of sensor -->
       <direction> 270.0 </direction>
   </timeSeries>
   

2. Exponential (power law) velocity profile:

   <!-- Start of timestep informastion for a sensor -->
   <timeSeries>
       <!-- time of the current data input -->
       <timeStamp>2003-07-23T23:00:00</timeStamp>
       <!-- Site BL flag (1-log (default), 2-exp, 3-urban canopy, 4-data entry) -->
       <boundaryLayerFlag> 2 </boundaryLayerFlag>
       <!-- Site exponent (using z0 variable) -->
       <siteZ0> 0.1 </siteZ0>
       <!-- Reciprocal Monin-Obukhov Length (1/m) -->
       <reciprocal> 0.0 </reciprocal>
       <!-- Height of the sensor -->
       <height> 20.0 </height>
       <!-- Measured speed at the sensor height -->
       <speed> 5.0 </speed>
       <!-- Wind direction of sensor -->
       <direction> 270.0 </direction>
   </timeSeries>
   

3. Urban canopy velocity profile:

   <!-- Start of timestep informastion for a sensor -->
   <timeSeries>
       <!-- time of the current data input -->
       <timeStamp>2003-07-23T23:00:00</timeStamp>
       <!-- Site BL flag (1-log (default), 2-exp, 3-urban canopy, 4-data entry) -->
       <boundaryLayerFlag> 3 </boundaryLayerFlag>
       <!-- Site z0 -->
       <siteZ0> 0.1 </siteZ0>
       <!-- Reciprocal Monin-Obukhov Length (1/m) -->
       <reciprocal> 0.0 </reciprocal>
       <!-- Height of the sensor -->
       <height> 20.0 </height>
       <!-- Measured speed at the sensor height -->
       <speed> 5.0 </speed>
       <!-- Wind direction of sensor -->
       <direction> 270.0 </direction>
       <!-- Height of the canopy -->
       <canopyHeight> 10.0 </canopyHeight>
       <!-- attenuation coefficient -->
       <attenuationCoefficient> 1.0 </attenuationCoefficient>
   </timeSeries>
   

4. Data entry of the profile from an experimental tower with multiple sensors or from a numerical mesoscale weather prediction model like WRF [13]:

   <!-- Start of timestep information for a sensor -->
   <timeSeries>
       <!-- time of the current data input -->
       <timeStamp>2003-07-23T23:00:00</timeStamp>
       <!-- Site BL flag (1-log, 2-exp, 3-urban canopy, 4-data entry) -->
       <boundaryLayerFlag> 4 </boundaryLayerFlag>
       <!-- Site z0 -->
       <siteZ0> 0.1 </siteZ0>
       <!-- Reciprocal Monin-Obukhov Length (1/m) -->
       <reciprocal> 0.0 </reciprocal>
       <!-- Height of the sensor -->
       <height> 30.7015 </height>
       <height> 74.4169 </height>
       <height> 144.644 </height>
       <height> 197.455 </height>
       <height> 268.468 </height>
       <!-- Measured speed at the sensor height -->
       <speed> 2.56922 </speed>
       <speed> 2.55532 </speed>
       <speed> 2.33319 </speed>
       <speed> 2.16058 </speed>
       <speed> 1.98843 </speed>
       <!-- Wind direction of sensor -->
       <direction> 323.283 </direction>
       <direction> 327.377 </direction>
       <direction> 332.676 </direction>
       <direction> 337.649 </direction>
       <direction> 344.273 </direction>
   </timeSeries>
   

Building Parameters (buildingsParams)

QES-Winds only conserves mass and no momentum equation is solved. As a result, the solution is a potential-flow solution (no shear effects). In order to add shear effects to our solution, empirical parameterizations are needed. These parameterizations are designed using results of experiments and computational simulations (e.g. [1, 6]). Buildings are the most important elements in urban areas. There are several parameterizations developed for different areas around the building. This section covers available parameterizations in QES-Winds along with their effects on the wind field.

<buildingsParams>
    <!-- Address to shapefile location-->
    <SHPFile>SaltLakeCity/slc_cut.shp</SHPFile>
    <!-- Name of building layer in shapefile-->
    <SHPBuildingLayer>slc_cut</SHPBuildingLayer>
    <!-- Name of building height field in shapefile -->
    <SHPHeightField>MEANHEIGHT</SHPHeightField>
    <!-- Height factor multiplied by the building height in the shapefile (default = 1.0)-->
    <heightFactor> 1.0 </heightFactor>

    <wallRoughness>0.01</wallRoughness>

    <!-- Upwind cavity flag (0-none, 1-Rockle, 2-MVP (default), 3-HMVP) -->
    <upwindCavityFlag> 2 </upwindCavityFlag>
    <!-- Wake flag (0-none, 1-Rockle, 2-Modified Rockle (default), 3-Area Scaled) -->
    <wakeFlag> 2 </wakeFlag>
    <!-- Street canyon flag (0-none, 1-Roeckle w/ Fackrel (default)) -->
    <streetCanyonFlag> 1 </streetCanyonFlag>
    <!-- Rooftop flag (0-none, 1-log profile (default), 2-vortex) -->
    <rooftopFlag> 1 </rooftopFlag>
    <!-- Sidewall flag (0-off, 1-on (default)) -->
    <sidewallFlag> 1 </sidewallFlag>
    <!--Street intersection flag (0-off (default), 1-on) -->
    <streetIntersectionFlag> 0 </streetIntersectionFlag>
    <!-- High-rise flag (0-off (default), 1-on) -->
    <highRiseFlag> 0 </highRiseFlag>
</buildingsParams>

Automated City Building

A new shapefile reader function has been added to QES-Winds, which provides the capacity to load the ESRI shapefiles using GDAL (Geospatial Data Abstraction Library) libraries. After the building footprints and heights are loaded from ESRI shapefiles, QES-Winds creates polygon buildings and applies appropriate parameterization to them. Figure below shows an example ESRI shapefile can be read into QES-Winds, Central Business District (CBD) of Oklahoma City shapefile, subject to JU2003 experimental campaign [14], plotted using the freely available software QGIS (https://qgis.org).

Central Business District (CBD) of Oklahoma City shapefile, subject to JU2003 experimental campaign [14], plotted using the freely available software QGIS.

The cell type contour for the Oklahoma City test case in a horizontal plane at \(z=3\) m is shown in Figure below. This plot indicates the ability of QES-Winds to read in and process ESRI shapefiles. The cell type value \(0\) (blue) represents the building cells while value \(1\) (red) indicates the air cells.

Cell type contour for the Oklahoma City test case in a horizontal plane at \(z=3\) m. The cell type value \(0\) (blue) represents the building cells while value \(1\) (red) indicates the air cells.

The user can define the address to the shapefile using <SHP> variable as well as the name of the shapefile using the <SHPBuildingLayer> and the correlation factor between the height field of the shapefile and the actual height of the buildings using the <heightFactor> under <simulationParameters> part in the XML file:

<buildingsParams>
    ...
    <!-- Address to shapefile location-->
    <SHP>../data/GISFiles/OKCSmallDomain/OKCSmallDomainJU2003.shp</SHP>
    <SHPBuildingLayer>OKCSmallDomainJU2003</SHPBuildingLayer>
    <!-- Height factor multiplied by the building height in the shapefile (default = 1.0)-->
    <heightFactor> 1.0 </heightFactor>
    ...
<buildingsParams>

Import Building From XML

Instead of reading in a ESRI shapefile, the user can import building information manually through the XML file. This can be done by using the <buildings> section of the XML file. The only option available for now is the rectangular building. Information required for defining a rectangular building are height, base height, length, width, location of the closest corner to the origin of domain and building rotational angle. Following is an example of a rectangular building with \(40\) m as height, \(0\) m as base height, \(20\) m as length and width, closest corner to the origin located at \(90\) m in \(x\) and \(y\) directions, and \(0^{\circ}\) as rotation angle with respect to the North-South line. Also, \(0.01\) m is defined as the surface roughness for all the building walls.

<buildingsParams>
    ...
    <wallRoughness> 0.01 </wallRoughness>
    <rectangularBuilding>
        <height> 40.0 </height>
        <baseHeight> 0 </baseHeight>
        <xStart> 90.0 </xStart>
        <yStart> 90.0 </yStart>
        <length> 20.0 </length>
        <width> 20.0 </width>
        <buildingRotation> 0.0 </buildingRotation>
    </rectangularBuilding>
    ...
<buildingsParams>

Upwind Cavity

Upwind cavity as described in [15, 16, 17] is the parameterization representing upwind and stagnation effects of the building on the fluid flow. There are three options available for this type of parameterization in QES-Winds.

The first option based on the parameterization proposed by Röckle [18] and later Kaplan and Dinar [19]. They defined an ellipsoid to represent what they call is the displacement zone in front of the building. The length of the displacement zone, \(L_F\), is defined by:

\[\frac{L_F}{H}=\frac{2(W / H)}{1+0.8 W / H} \label{eq:lf}\]

The shape of the ellipsoid is estimated by:

\[\frac{X^{2}}{L_F^{2}\left(1-(Z / 0.6 H)^{2}\right)}+\frac{Y^{2}}{W^{2}}=1 \label{eq:upwind}\]

where \(L\), \(H\) and \(W\) are length, width and height of the building, receptively.Finally, the initial velocity components in the displacement zone are set to zero.

Part (a) of figures below show cell type contour to represent the area of effect of the Röckle upwind cavity parameterization in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively. The upwind parameterizations is applied to a rectangular building with the initial guess field is constructed using a single sensor with logarithmic profile. Parts (b) and (c) of figures below indicate velocity magnitude contour with overlaying velocity vectors of initial (part (b)) and final (part(c)) velocity fields in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively.

(a) Cell type contour to show the area of effect of the Röckle upwind cavity parameterization in a vertical plane at \(y=100\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a vertical plane at \(y=100\) m.

(a) Cell type contour to show the area of effect of the Röckle upwind cavity parameterization in a horizontal plane at \(z=5\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a horizontal plane at \(z=5\) m.

The second option is called the Modified Vortex Parameterization (MVP) and created by Bagal et al. [16]. In this parameterization, the length of the displacement zone, \(L_F\), is calculated by equation below. The MVP parameterization defines two ellipsoids instead of one: In the outer ellipsoid, velocities are reduced to \(40\%\) of their initial values while in the inner region, velocity components are set to zero [15]. Both ellipsoids are extended to \(0.6\) of the building height.

\[\frac{L_F}{H}=\frac{1.5(W / H)}{1+0.8 W / H} \label{eq:lf_MVP}\]

where \(L\), \(H\) and \(W\) are length, width and height of the building, receptively.

Part (a) of figures below show cell type contour to represent the area of effect of the MVP upwind cavity parameterization in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively. The upwind parameterizations is applied to a rectangular building with the initial guess field is constructed using a single sensor with logarithmic profile. Parts (b) and (c) of the figures below indicate velocity magnitude contour with overlaying velocity vectors of initial (part (b)) and final (part(c)) velocity fields in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively.

(a) Cell type contour to show the area of effect of the MVP upwind cavity parameterization in a vertical plane at \(y=100\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a vertical plane at \(y=100\) m.

(a) Cell type contour to show the area of effect of the MVP upwind cavity parameterization in a horizontal plane at \(z=5\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a horizontal plane at \(z=5\) m.

The third option is called the high-rise MVP algorithm (HMVP) and is designed to address the shortcomings of the previous models when it comes to tall buildings [15]. The length of the displacement zone \(L_F\) is calculated using the equation presented above. The HMVP algorithm creates two ellipsoids with the difference that the inner region only extends to \(60\%\) of the minimum of building height and building width. In addition, the algorithm linearly reduces the velocities in the outer region from their upwind values at the outer surface to \(40\%\) of the initial values on the inner region.

Part (a) of figures below show cell type contour to represent the area of effect of the HMVP upwind cavity parameterization in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively. The upwind parameterization is applied to a rectangular building with the initial guess field is constructed using a single sensor with logarithmic profile. Parts (b) and (c) of the figure below indicate velocity magnitude contour with overlaying velocity vectors of initial (part (b)) and final (part(c)) velocity fields in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively.

(a) Cell type contour to show the area of effect for the HMVP upwind cavity parameterization in a vertical plane at \(y=100\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a vertical plane at \(y=100\) m.

(a) Cell type contour to show the area of effect of the HMVP upwind cavity parameterization in a horizontal plane at \(z=5\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a horizontal plane at \(z=5\) m.

In order to choose between these three upwind models, the user needs to change the value of “upwindCavityFlag” in the XML file.

<buildingsParams>
    ...
    <!-- Upwind cavity flag (0-none, 1-Rockle, 2-MVP (default), 3-HMVP) -->
    <upwindCavityFlag> 2 </upwindCavityFlag>
    ...
</buildingsParams>

Leeside Cavity and Far-Wake

The far-wake and cavity parameterization described in [20, 21] are a significant part of the building parameterizations. The one available in QES-Winds is based on the parameterization proposed by Röckle [18] and later Kaplan and Dinar [19]. The Röckle parameterization defines two ellipsoids to represent the shape of the reversed flow cavity and the far-wake region. The reversed flow cavity extends to the along-wind cavity length (\(L_R\)), which is calculated as:

\[\frac{L_{R}}{H}=\frac{1.8 \frac{W}{H}}{\left(\frac{L}{H}\right)^{0.3}\left(1+0.24 \frac{W}{H}\right)}, \label{eq:Lr}\]

and wake is assumed to be approximately \(3\) cavity lengths long (i.e., \(3L_R\)). After calculating \(L_R\), the cavity length, \(d\) in the stream-wise direction was defined by an ellipsoid shape using:

\[d=L_{R} \sqrt{\left(1-\left(\frac{z}{H}\right)^{2}\right)\left(1-\left(\frac{y}{W}\right)^{2}\right)}-\frac{L}{2}. \label{eq:d}\]

Finally, the velocity in the reversed cavity zone is defined using:

\[\frac{u(x, y, z)}{U(H)}=-\left(1-\left(\frac{x}{d}\right)^{2}\right) \label{eq:cavity}\]

and in the wake region, the velocity field is estimated by:

\[\frac{u(x, y, z)}{U(H)}=\left(1-\left(\frac{d}{x}\right)^{1.5}\right). \label{eq:wake}\]

where \(L\), \(H\) and \(W\) are length, width and height of the building, receptively. \(u(x,y,z)\) is the velocity at point \((x,y,z)\), \(U(H)\) is the reference velocity at height of the building and \(x\) is the distance from the building in the stream-wise direction.

Part (a) of the figure below show cell type contour to represent the area of effect of the Röckle wake parameterization in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively. The wake parameterization is applied to a rectangular building with the initial guess field is constructed using a single sensor with logarithmic profile. Parts (b) and (c) of the figures below indicate velocity magnitude contour with overlaying velocity vectors of initial (part (b)) and final (part(c)) velocity fields in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively.

(a) Cell type contour to show the area of effect of the Röckle wake parameterization in a vertical plane at \(y=100\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a vertical plane at \(y=100\) m.

(a) Cell type contour to show the area of effect of the Röckle wake parameterization in a horizontal plane at \(z=5\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a horizontal plane at \(z=5\) m.

In order to turn on the wake model, the user needs to change the value of “wakeFlag” in the XML file.

<buildingsParams>
    ...
    <!-- Wake flag (0-none, 1-Rockle, 2-Modified Rockle (default), 3-Area Scaled) -->
    <wakeFlag> 2 </wakeFlag>
    ...
<buildingsParams>

Street Canyon

The street canyon parameterization detailed in [6] represents the effects of two buildings in close vicinity to each other, on the fluid flow. Röckle [18] introduced velocity parameterizations for the stream-wise components as:

\[\frac{u(x, y, z)}{U(H)}=-\frac{x_{\mathrm{can}}}{(0.5 S)}\left(\frac{S-x_{\mathrm{can}}}{0.5 S}\right) \label{eq:u_can}\]

and the vertical component as

\[\frac{w(x, y, z)}{U(H)}=-\left|\frac{1}{2}\left(1-\frac{x_{\text {can }}}{0.5 S}\right)\right|\left(1-\frac{S-x_{\text {can }}}{0.5 S}\right). \label{eq:w_can}\]

where \(S\) is the spacing between two buildings and \(x_{can}\) is the distance from the backwall of the upwind building.

In order to identify the criteria to determine the existence of a street canyon, Singh et al. [6] utilized the cavity length, \(L_R\), for the upwind building. If \(S < L_R\), the street canyon parameterization is applied, otherwise, the upwind building is considered as an isolated building.

Part (a) of the figures below show cell type contour to represent the area of effect of the street canyon parameterization in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively. The street canyon parameterization is applied to an area between two rectangular buildings. The upwind building is same as the one defined previously. The downwind building is a rectangular building with \(20\) m as height, \(0\) m as base height, \(20\) m as length and width, closest corner to the origin located at \(90\) m in \(x\) and \(120\) m in \(y\) directions, and \(0^{\circ}\) as rotation angle with respect to the North-South line. The initial guess field is constructed using a single sensor with logarithmic profile. Parts (b) and (c) of the figures below indicate velocity magnitude contour with overlaying velocity vectors of initial (part (b)) and final (part(c)) velocity fields in a vertical plane at \(y=100\) m and a horizontal plane at \(z=5\) m, respectively.

(a) Cell type contour to show the area of effect of the street canyon parameterization in a vertical plane at \(y=100\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a vertical plane at \(y=100\) m.

(a) Cell type contour to show the area of effect of the street canyon parameterization in a horizontal plane at \(z=5\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a horizontal plane at \(z=5\) m.

To turn on the street canyon parameterization, the user needs to change the value of “streetCanyonFlag” in the XML file.

<buildingsParams>
    ...
    <!-- Street canyon flag (0-none, 1-Roeckle w/ Fackrel (default)) -->
    <streetCanyonFlag> 1 </streetCanyonFlag>
    ...
<buildingsParams>

Rooftop Recirculation

The rooftop parameterization described in [22, 23], captures the separation of the flow from the leading edge of the building. It first checks if the incident flow is in \(\pm15^{\circ}\) of perpendicular to the front face. The parameterization then creates an ellipsoidal region above the building with height of \(H_c\) (height of the vortex) and length of \(L_c\) (length of the vortex). It applies a logarithmic profile in the whole vortex area and finally, reverses the velocity in region \(1\). Region \(1\) is an ellipsoidal zone with the same length as the vortex and half of the height.

\[R=B_{\mathrm{s}}^{2 / 3} B_{l}^{1 / 3}\]

\[L_{\mathrm{c}}=0.9 R \label{eq:Lc}\]

\[H_{\mathrm{c}}=0.22 R \label{eq:Hc}\]

where \(B_s\) is the smaller of the height (\(H\)) and the effective width (\(W_{eff}\)) of the building, \(B_l\) is the larger of \(H\) and \(W_{eff}\) , \(R\) is the vortex size scaling factor.

Part (a) of the figure below show cell type contour to represent the area of effect of the rooftop parameterization in a vertical plane at \(y=100\) m. The rooftop parameterization is applied to a rectangular building with \(40\) m as height, \(0\) m as base height, \(40\) m as length and width, closest corner to the origin located at \(90\) m in \(x\) and \(y\) directions, and \(0^{\circ}\) as rotation angle with respect to the North-South line. The initial guess field is constructed using a single sensor with logarithmic profile. Parts (b) and (c) of the figure below indicate velocity magnitude contour with overlaying velocity vectors of initial (part (b)) and final (part(c)) velocity fields in a vertical plane at \(y=100\) m.

(a) Cell type contour to show the area of effect of the rooftop parameterization in a vertical plane at \(y=100\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a vertical plane at \(y=100\) m.

To turn the parameterization on, the user needs to change the value of “rooftopFlag” in the XML file.

<buildingsParams>
    ...
    <!-- Rooftop flag (0-none, 1-log profile (default)) -->
    <rooftopFlag> 1 </rooftopFlag>
    ...
<buildingsParams>

Sidewall Recirculation

The sidewall parameterization is designed to represent the effects of the edge of the building on the upwind field [4]. It first checks if a face has an outward normal vector nominally (\(\pm 10^{\circ}\)) perpendicular to the local wind vector. The important parameters controlling the sidewall vortex strength and geometry are:

\[R=B_{\mathrm{s}}^{2 / 3} B_{l}^{1 / 3}\]

\[L_{\mathrm{c}}=0.9 R\]

\[W_{\mathrm{c}}=0.22 R\]

where \(B_s\) is the smaller of the height (\(H\)) and the effective width (\(W_{eff}\)) of the building, \(B_l\) is the larger of \(H\) and \(W_{eff}\) , \(R\) is the vortex size scaling factor, \(L_c\) is the downwind length of the half-ellipse that defines the vortex recirculation region, and \(W_c\) is the lateral width of the elliptical recirculation region. Within the recirculation zone, the velocity is reversed and scaled linearly from the reference wind speed near the wall to zero at the edge of the ellipse.

Part (a) of the figure below show cell type contour to represent the area of effect of the sidewall parameterization in a horizontal plane at \(z=5\) m. The rooftop parameterization is applied to a rectangular building with the initial guess field is constructed using a single sensor with logarithmic profile. Parts (b) and (c) of the figure below indicate velocity magnitude contour with overlaying velocity vectors of initial (part (b)) and final (part(c)) velocity fields in a horizontal plane at \(z=5\) m.

(a) Cell type contour to show the area of effect of the sidewall parameterization in a horizontal plane at \(z=5\) m. Velocity magnitude contour with overlaying velocity vectors of (b) initial velocity field and (c) final velocity field, in a horizontal plane at \(z=5\) m.

In order to turn the algorithm on, the user needs to change the value of “sidewallFlag” in the XML file.

<buildingsParams>
    ...
    <!-- Sidewall flag (0-off, 1-on (default)) -->
    <sidewallFlag> 1 </sidewallFlag>
    ...
<buildingsParams>

Street Intersection

This parameterization is in developpement

In order to turn the parameterization on, the user needs to change the value of “streetIntersectionFlag” in the XML file.

<buildingsParams>
    ...
    <!--Street intersection flag (0-off (default), 1-on) -->
    <streetIntersectionFlag> 0 </streetIntersectionFlag>
    ...
<buildingsParams>

High-rise Parameterization

This parameterization is in developpement

In order to turn the parameterization on, the user needs to change the value of “highRiseFlag” in the XML file.

<buildingsParams>
    ...
    <!-- High-rise flag (0-off (default), 1-on) -->
    <highRiseFlag> 0 </highRiseFlag>
    ...
<buildingsParams>

Vegetation Parameters (vegetationParams)

Row-organized canopy (ROC) model

The ROC model adjusts the mean wind field to account for drag in sparse, structured row crops (e.g., grape vineyards, carrots, and some orchards). It is comprised of several parameterizations that alter the flow in specific regions around each row in the ROC. These zones are pictured in Figure 1.1. The ROC model is described in detail in [24].

Schematic of ROC model. (a) Side view of ROC with a sparse understory, (b) Side view of ROC without understory space. Green patches indicate vegetation (or fence) where the bleed flow parameterization is applied. Red patches indicate the region where the wake parameterization is applied. The entry region is shown in light red, and the equilibrium and recovery regions are shown in dark red. Black lines behind the last row in (a) demonstrate the shape of the sheltering profile (sparse understory) multiplied with a logarithmic inlet velocity profile. The black lines in (b) demonstrate the no-understory version. Blue patches indicate the high-pressure upwind deficit (UD) zones. Dotted lines show the boundaries of the shear mixing zones. (c) Top-down view of the object-aligned coordinate system for the ROC model. The coordinate \(x'\) is the streamwise distance downwind of each row, and is aligned with the reference wind direction and \(x\) is the row-orthogonal distance from a row.

A wake model forms the basis of the ROC parameterizations. It applies a wake from each row of vegetation separately, and the superposition of all wakes forms the overall canopy flow. To approximate the effect of drag from the vegetation, the initial velocities behind each row are multiplied by a sheltering coefficient \(\alpha(x,y,z) \leq 1\). An algebraic expression is prescribed for \(\alpha\) that mimics the shape of mean velocity profiles found behind windbreaks [25]. At any point in the canopy, \(\alpha\) contains contributions from the nearest upwind row (“local” sheltering) and all other rows upwind of it (“non-local” sheltering). The sheltering profile from a single row is:

\[\label{wakeAlpha} \alpha_r(x',z) = \frac{1-\alpha_{o}}{2}\tanh\left(1.5\frac{z-z_{szo}^{top}}{\lambda^{top}}\right) + \frac{1+\alpha_{o}}{2}\]

for \(z \geq z_{mid}\), where \(z_{mid}\) is the vertical midpoint of the canopy, and

\[\alpha_r(x',z) = \frac{1-\alpha_{o}}{2}\tanh\left(1.5\frac{z_{szo}^{bot}-z}{\lambda^{bot}}\right) + \frac{1+\alpha_{o}}{2}\]

for \(z < z_{mid}\). The superscripts \(top\) and \(bot\) refer to separate shear zones at the canopy top and understory, respectively, and \(\lambda\) is the thickness of a given shear mixing layer. The shear zone thicknesses increase linearly with downwind distance from each row (the coordinate \(x'\) denotes the distance from the leeward side of a row along a vector aligned with the reference wind direction). Growth of \(\lambda\) causes the vertical gradient of \(\alpha(x',z)\) to decrease with \(x'\), representing wake recovery behavior. The growth rate of \(\lambda\) (i.e., the recovery speed) increases with background turbulence (quantified by \(\sigma_w\)), which causes the mixing layer to diffuse vertically. It also increases with vegetation density; with denser vegetation, the velocity gradient driving the growth of the mixing layer increases. The vertical center (i.e., origin) of a given shear mixing zone is denoted by \(z_{szo}\). This origin gradually drops to the ground as the wake recovers. This movement is parameterized with experimental data from [26]. The quantity \(\alpha_o\) is the aerodynamic porosity of a single row and is obtained from the leaf area density (LAD) of the row. The aerodynamic porosity forms the primary boundary condition of the wake parameterization.

The total sheltering behind any given row in the ROC is not simply the product of all \(\alpha\) terms from all upwind rows; the velocity would approach zero for large row numbers. Instead, the ROC is partitioned into two regions: an entry region at the windward edge of the ROC, where momentum is lost from row to row, and an equilibrium region, where velocities are periodically oscillating. In the equilibrium zone, momentum is maintained from row to row as losses to canopy drag are balanced by turbulent momentum flux from above. In the model, these zones are delineated through different constructions of the total sheltering. In the entry region, the total sheltering behind any row is the product of the \(\alpha\) terms for all rows upwind of it. In the equilibrium region, the sheltering is the product of \(\alpha\) from the nearest row (the “local” row) and the rows in the entry region only (“non-local” rows), which caps the maximum net momentum loss within the canopy. Calculation of the number of rows in the entry region is explained in [24].

Upwind of each row, there is a region of high pressure where momentum is lost to form drag. In the model, a parabolic region attached to the windward side of each row is defined in which velocities are reduced and recirculation may occur. Within that region, a velocity deficit is calculated based on a scaling argument that the vegetative drag term in the mean momentum equation is of similar magnitude as the advection term:

\[\label{UDscaling} C_d \cdot LAD_{avg} \cdot U^2 \sim U \frac{\partial U}{\partial x}\]

The velocity deficit, which originates from the \(\partial U\), is subtracted from the existing wake velocities throughout the parabolic region.

ROCs are added to a QES domain as a polygonal block, rather than as individual rows. Polygon vertices are provided by the user and rows are automatically generated within the block according to the user-specified row geometry parameters (spacing, thickness, height, understory height). The XML for adding ROCs is as follows:

<vegetationParams>
  <vegetationParams>
        <num_canopies>1</num_canopies>
        <ROC>
            <opticalPorosity>0.197</opticalPorosity>
            <!-- Height of canopy -->
            <height>2.16</height>
            <!-- Height at which vegetation begins. Zero if the canopy extends to the ground. -->
            <understoryHeight>0.5</understoryHeight>
            <!-- Front-to-front distance between rows -->
            <rowSpacing>2.5</rowSpacing>
            <!-- Thickness of a single row. (still required for thin-fence parameterization) -->
            <rowWidth>0.5</rowWidth>
            <!-- Row angle: orientation of rows, in compass degrees -->
            <rowAngle>0</rowAngle>
            <xVertex>5</xVertex>
            <xVertex>45</xVertex>
            <xVertex>45</xVertex>
            <xVertex>5</xVertex>
            <yVertex>5</yVertex>
            <yVertex>5</yVertex>
            <yVertex>55</yVertex>
            <yVertex>55</yVertex>
            <!-- Thin fence flag: 1 for thin fence (aero. porosity = optical porosity), 2 for vegetative row (aero. por. from LAD) -->
            <thinFence>0</thinFence>
            <!-- Upstream sig_w used in wake recovery rate calculation -->
            <upstreamSigmaW>0.5187</upstreamSigmaW>
            <!-- Upstream ustar used in entry length calculation -->
            <upstreamUstar>0.5067</upstreamUstar>
            <!-- Displacement height used in row-parallel component (empirical) -->
            <displacementHeightParallel>1.42</displacementHeightParallel>
            <!-- Average LAD: vertically averaged LAD for a single row, used in upwind displacement zone drag calculation -->
            <LAD_avg>4.6418</LAD_avg>
            <!-- Effective LAD: spatially averaged LAD (as in Bailey & Stoll, 2013), a bulk drag parameter used in entry length calculation -->
            <LAD_eff>0.71736</LAD_eff>
            <!-- TKE above the canopy, used in the k-l turbulence model -->
            <tkeMax>1.9746</tkeMax>
        </ROC>
    </vegetationParams>