Effect of different source terms in atmospheric boundary modelling over the complex terrain site of Perdigao

The assessment of wind conditions in complex terrain requires the use of Computational Fluid Dynamics (CFD) simulations incorporating an accurate parameterization of forest canopy effects and variable thermal stability effects. This study aims to investigate how incorporating the presence of trees can improve the flow predictions. A three-dimensional steady Reynolds-averaged Navier-Stokes (RANS) equations model is set up using OpenFOAM to simulate the flow over a complex terrain site comprising two parallel ridges located near Perdigão, Portugal. A 7.5 km × 7.5 km terrain of the Perdigao site is 5 constructed from elevation data centered around a 100 m met-mast located on the northeast ridge. A 30-min averaged stationary period corresponding to near-neutral conditions on a single met-mast tower is simulated. The impact of incorporating different source terms is studied such as forest canopy, Coriolis forces as well as also buoyancy forces. The prediction capability of the models is analyzed for different groups of towers on the South-West ridge, inside the valley and on the North-East ridge based on the flow topology. The inclusion of a canopy model is shown to improve predictions close to the ground for most of the 10 towers, while reducing prediction accuracy on top of the ridges, illustrating the need to represent terrain heterogeneity.

flows over complex terrains to account for such phenomena and provide more reliable wind resource predictions (Blocken, 2014). However, improving numerical modeling tools for complex terrains demands accounting for various microscale phe-25 nomena such as flow re-circulation and the forest canopy effect. Such microscale flow features have a significant impact on local wind resource assessment and wind turbine loading.
Inital studies in literature accounting for effects of forest canopy were performed over forested flat terrains. Brower (2012) showed that the presence of a forest canopy increases the modelling uncertainty by up to a factor of 5 irrespective of the chosen modelling approach. Finnigan (2000) parameterized the effect of foliage and forest canopy by accounting for the drag force in 30 the momentum equation. The effect of canopy on turbulence was further taken into account by Sogachev and Panfyorov (2006) and Sogachev (2009) by adding additional source terms for turbulence kinetic energy and turbulence dissipation rate. Desmond et al. (2017) modelled both forest canopy and buoyancy effects by modelling them using source and sink terms and showed that the thermal stratification plays an important role in determining the flow over canopies. The canopy model is typically implemented by specifying the tree height and leaf area density (LAD), which represents the area of the leaf and branches of trees. 35 These parameters are typically defined from field measurements or LiDAR point clouds generated from field measurement campaigns and aerial surveys (Queck et al., 2012).
Validation of numerical models for complex terrains requires field measurements. One of the earliest field campaigns was at Askervin hill (Salmon et al., 1988), which is a smooth isolated hill. Similarly, the Bolund hill campaign was another campaign over an isolated hill in Denmark (Bechmann et al., 2009) with data under conditions which could be classified as neutral. 40 However no canopy is seen on these sites. The Alaiz field campaign (Rodrigo et al., 2021) was performed over a large scale homogeneous mountain valley topography in Spain with forested regions. Chavez et al. (2014) showed that using a canopy model improved the accuracy of velocity and turbulence predictions for simulations over Alaiz.
More recently, a field experiment campaign has been carried out over the complex terrain terrain site at Perdigão, Portugal ( (Fernando et al., 2018)) which is a double ridge with forested canopy experiencing different thermal stratification condi-45 tions. Significant field data is available for further studies and is further elaborated in Section 2. Several numerical studies have been performed for validation with field data. Laginha Palma et al. (2020)  This was chosen due to practical reasons related to the model setup and is not representative of actual tree height. There is a 50 gap in existing literature on microscale simulations using a canopy model on a heterogeneous forested complex terrain with a rich dataset such as Perdigão. This study aims to address that research gap by investigation of the effect of specific microscale effects related to the modeling of trees and buoyancy effects at the complex terrain test-case of Perdigão at several met-mast locations and to verify if the forest parameterization helps improve predictions using a RANS simulation approach.

55
This paper aims to evaluate the impact of different physical source terms and turbulence models when performing CFD simulations in complex terrain and compare simulation prediction with field measurement profiles at the different groups of towers      Easting ( (Kaimal and Finnigan, 1994) is calculated using Eqn 1. However, the flow conditions are non-neutral at other met mast 75 locations.
where θ is the potential temperature, g is the acceleration due to gravity, U, V are the wind velocity components at the reference height (m). The same period was simulated by Laginha Palma et al. (2020) based on the stationary periods predicted by Carvalho (2019).

80
The terrain is a 7.5 km × 7.5 km square centralized around a 100 m met-mast located on the south-west ridge.  A cylindrical computation domain was developed, which provides the flexibility to simulate wind from any wind direction. This approach has also been successfully applied by the authors in a previous study pertaining to urban areas (Hågbo 85 et al. 2020). A smoothing region from complex to flat terrain was applied towards the outer boundaries, with a minimum radial distance of 15 ∆h t . Several best practice guidelines have been formulated for grid generation for simulating complex terrains, such as by Sørensen et al. (2012) and Laginha Palma et al. (2020) and have been closely followed. The height of the domain is set to ten times the difference in the elevation height of the terrain, ∆h t , as recommended by Sørensen et al. (2012) when simulating wind flow over complex terrains. It consists of 12.7 million cells and is produced with 90 terrainBlockMesher (Schmidt et al., 2012), capable of generating structured meshes over complex terrain exclusively consisting of hexahedra cells. Terrestrial data was obtained from the Shuttle Radar Topography Mission (SRTM) database of the Perdigão field experiment (Fernando et al., 2018). The horizontal mesh resolution over the terrain was set to 33 m. It satisfies the minimum resolution of 40 m recommended specifically for the Perdigão site by Laginha Palma et al. ( 2020). A uniform stretching is applied in the vertical direction. A mesh sensitivity study was performed to ensure that the grid resolution  the finite volume method. Second-order discretization schemes were used for the spatial discretization. The initial iterative convergence criteria were that the scaled residuals should drop four orders of magnitude for all flow variables as per the BPGs. Two steady-state solvers for turbulent flow of incompressible fluids have been used: simpleFoam (SF) and buoyantBoussinesqSimpleFoam (BBSF). All thermal effects are neglected in the simulations using the SF solver, such that the atmospheric stability of the simulated atmosphere is always neutral. In these simulations, the air density is assumed 105 constant, and the gravitational force is neglected. The effect of buoyancy forces is included in the simulations using the BBSF solver. In addition to solving the continuity equation and the momentum equation, which is solved using SF, the energy equation is solved. It allows for the modeling of non-neutral atmospheric conditions. The solver uses the Boussinesq approximation, which ignores density differences except where they appear in terms multiplied by the gravitational acceleration.

Inlet profiles and Boundary conditions 110
Two sets of fixed inlet profiles have been used, both forming a homogeneous atmospheric boundary layer (ABL). either idealized with a logarithmic velocity profile or fully developed profiles obtained from a precursor simulation. The precursor model is further described by Alletto et al. (2018) and Koblitz (2013). The idealized profiles provide the turbulence quantities and a log-law type ground-normal wind velocity based on the generalization and modification of the well-used set of equations from Richards and Hoaxey (1993). This modification has been 115 implemented by Yang et al. (2009) and uses a more mathematically consistent formulation allowing the turbulent kinetic energy to vary with height. Also, using a log-low wind velocity profile and the associated turbulent inflow conditions are considered only to be valid in the atmospheric surface layer, which can be roughly estimated as 10% of the atmospheric boundary layer (Temel et al. 2018).
An alternative to using the idealized profiles is using fully developed profiles obtained from a one-dimensional precursor to produce inlet profiles valid for various conditions, including non-neutral atmospheric conditions. Adjusting the parameters of these source terms is an efficient way to obtain the desired inflow. An overview of the boundary conditions (BCs) used is presented in Table 1.  applied ( Hargreaves and Wright, 2007). The roughness was set to z 0 = 0.02 m, based on average values of a roughness map provided at the database of the Perdigão field experiment (Fernando et al., 2018).

Source terms
Four different source terms have been added to some of the simulation cases: Tree canopy as a porous medium, Coriolis force, Pressure gradient force, Buoyancy forces, maximum mixing length scale limiter for the mixing length turbulence model.

140
Some simulation cases include the modeling of a tree canopy through a source term. The canopy is modeled as a porous medium based on the work of Costa ( 2007). A drag term is added to the momentum equation and additional source terms for the production and dissipation of turbulence. The porous media is implemented through the power law porosity model. The parameters for the canopy modeling are the leaf area density, which is set to 0.14 m −1 , and the drag coefficient of the trees, where LAD is the Leaf area density, C d is the plant canopy drag coefficient. The magnitude and direction of the Coriolis force source term is calculated based on Earths' rotation and the Coriolis frequency parameter corresponding to the latitude of Perdigão shown in Eq. 3.
where Ω is the rotational period and ρ is the air density.
A pressure gradient force is also included in the simulation cases using the BBSF solver. In setups using this solver, the wind 155 is driven by the balance of the lateral pressure gradient force and the Coriolis force. This balance is called geostrophic balance, and the resulting wind is a geostrophic wind, which in the Northern Hemisphere goes counterclockwise around low-pressure systems.
The inclusion of buoyancy forces is enabled by adding a heat/cooling flux on the terrain patch and adjusting the maximum mixing length scale, L max , through the use of the source terms for turbulent kinetic energy and dissipation shown in Eqn. 4.

160
These buoyancy forces are required in modeling non-neutral atmospheric effects and are only possible using the BBSF solver of the two solvers used .
where B is the buoyancy production term, β B is the thermal expansion coefficient, α to is the kinematic turbulent thermal conductivity, k is the turbulent kinetic energy, is the turbulent kinetic energy dissipation rate, C 3 is a model constant. Terms 165 with the subscript denote values for previous iteration.
The level of turbulent kinetic energy is controlled in the domain by the use of a parameter called the maximum turbulent mixing length scale, used is to limit the production and dissipation terms of the turbulence model. The maximum global turbulent mixing length scale for a neutral ABL as is shown in Eq. 5 proposed by (Blackadar, 1962) and also used by (Koblitz, 2013).
where G is the geostrophic wind and f c is the Coriolis frequency parameter.

Simulation cases
The setup for all the simulations was identical except for certain parameters. The solver used is either simpleFoam (SF) or buoyantBoussinesqSimpleFoam (BBSF). The inlet profiles, either idealized with a logarithmic velocity profile repre- model was tuned using experimental data and LES (Large Eddy Simulation) for a homogeneous forest (Costa 2007). With regards to the boundary conditions, they are overall the same but had to be adjusted to account for the different solvers and inlet profiles used. An overview of all the simulation cases is provided in Table 2.

Results
The  Table 3. A hit rate metric for a given model is defined as the number of predictions predictions within one standard deviation over each height of the field measurement. The relative 190 error for the turbulent kinetic energy and wind direction is shown in Tables 4, 5 respectively. The relative error is defined as the percentage of difference between the simulated and field measurement, which is divided by field measurement value.

Model prediction on SW ridge
The towers on the SW ridge as shown in Fig. 2 are the Tower 20, 34 and 37. This is a region of flow acceleration. A good match is obtained in between the measured and computed velocity, turbulent kinetic energy and wind direction profiles for the 195 calibration Tower 20 as seen in Fig. 6. For Tower 34 shown in Fig 7 and for Tower 37 shown in Fig 8, the model including the canopy provides the best match. In models that do not account for a canopy, a speed up is seen close to the ground with large RMSE errors as seen in Table 3. The canopy model on the other hand predicts the shape of the profile correctly, but over predicts the turbulent kinetic energy for all the towers on the SW ridge. The results suggest that the inclusion of the canopy is a good choice for prediction of velocity profiles, however the parameter choice needs to be better optimised to reduce the 200 levels of turbulence produced due by the canopy. In terms of the turbulence model, the k − model is seen to produce a slightly    The measured values of turbulent kinetic energy is high inside the valley, indicative of flow mixing and high turbulence, however all the models seem to strongly under predict the turbulent kinetic energy with large relative errors as seen in Table 4, except for Tower 22 which is the closest to the South-West ridge. Almost all models seen to provide a good prediction of the velocity profiles and TKE at that location, however a sudden shift in wind direction, indicative of a flow separation and 210 re-circulation is seen in Fig 11. The buoyancy model provides a good hit rate for the wind direction inside the valley as shown in Table 5. On comparison of the turbulence the turbulence model, the k − model performs better than the k − ω model in prediction of wind velocity. Fig. 13 shows the velocity contour for the predicted re-circulation zone inside the valley using the canopy and buoyancy models. The model with the canopy produces a larger re-circulation zone compared to the other models. Other models over-predict the velocity profile close to the ground. However, the canopy parameters need to be tuned as the surface heterogeneity is not taken into account as the prediction accuracy varies at the different locations along the ridge. 235 b) For the towers inside the valley: The valley is a zone of flow re-circulation. and comprises lower velocities and higher variability which remains challenging for prediction models. The BBSF model is shown to provide the best prediction of wind direction. Large relative errors in wind speed and turbulent kinetic energy profiles are seen for most of the models. c) For the towers on the North-East ridge: The region at the North-East ridge is a zone of flow acceleration downstream 240 of the re-circulation zone from the valley. Significant turbulence is seen close to ground in the field measurements, which is under-predicted by most of the models. The model accounting for buoyancy forces provides a good prediction on Tower 29 for the wind speed and direction while the canopy model strongly under predicts the wind speed. Local terrain features need to be resolved such as the mountain gap near Tower 10.

245
Finally, the choice of the best turbulence model is found to be inconclusive in terms of overall prediction capability for different parts of the terrain. In the future, the surface heterogeneity of canopy is to be modelled based on the surface roughness map at Perdigão. The inclusion of the canopy model in buoyantBoussinesq solver could help evaluate the benefit of incorporating both canopy and buoyancy effects together. Furthermore, simulations shall be performed for periods of non-neutral flow conditions, such as stable or unstable thermal stratification conditions, where the buoyancy effects play an important role.