The digital terrain model in the computational modelling of the flow over the Perdigão site: the appropriate grid size

The digital terrain model (DTM), the representation of earth's surface at regularly spaced intervals, is the first input in the computational modelling of atmospheric flows. The ability of computational meshes based on high- (2 m; airborne laser scanning, ASL), medium- (10 m; military maps, Mil) and low-resolution (30 m; Shuttle Radar Topography Mission, SRTM) DTMs to replicate the Perdigão experiment site was appraised in two ways: by their ability to replicate the two main terrain attributes, elevation and slope, and by their effect on the wind flow computational results. The effect on the flow modelling was evaluated by comparing the wind speed, wind direction and turbulent kinetic energy using VENTOS^®/2 at three locations, representative of the wind flow in the region. It was found that the SRTM was not an accurate representation of the Perdigão site. A 40 m mesh based on the highest-resolution data yielded an elevation error of less than 1.4 m and an RMSE of less than 2.5 m at five reference points compared to 5.0 m in the case of military maps and 7.6 m in the case of the SRTM. Mesh refinement beyond 40 m yielded no or insignificant changes on the flow field variables, wind speed, wind direction and turbulent kinetic energy. At least 40 m horizontal resolution – threshold resolution – based on topography available from aerial surveys is recommended in computational modelling of the flow over Perdigão.

1 Introduction

A digital terrain or digital elevation model (DTM or DEM) is a representation of the earth's surface elevation at regularly spaced horizontal intervals. Although the terrain model is the first input in computational modelling of atmospheric flows, its impact on flow results has not been a matter of concern because the spatial resolution of publicly available DTMs is higher than the size of the computational grid often used to resolve the terrain. However, before a fluid flow database (Mann et al., 2017) can be used as a reference in flow model appraisal and development, the impact of the terrain modelling must be assessed. For studies of the atmospheric flow over Perdigão the publicly available DTMs were considered not accurate enough (Mukherjee et al., 2013; Simpson et al., 2015), and an airborne laser scanning (ASL) campaign of the region was carried out in 2015, first to assist the design of the Perdigão campaigns in 2015 and 2017 (cf., Vasiljević et al., 2017; Fernando et al., 2019) and second to provide the high-resolution terrain data for computational flow modelling, on par with the resolution provided by the large amount of measuring equipment within a small region.

The Perdigão site is located in the municipality of Vila Velha de Ródão, in the centre of Portugal (4396621N, 608250E: ED50 UTM 29N or, in WGS84, 39^∘42^′38.5^′′ N 7^∘44^′18.5^′′ W). It is comprised of two parallel ridges at about 500 m elevation, 4 km in length, with a SE–NW orientation and distanced around 1.4 km from each other. The land is covered by a mixture of farming areas and patches of eucalyptus (Eucalyptus globulus) and pine trees (Pinus pinaster). The dominant winds are perpendicular to the ridges, assuring a largely two-dimensional flow.

The accuracy of a DTM depends on the data collection techniques, data sampling density and data post-processing, such as grid resolution and interpolation algorithms. In computational modelling of atmospheric flows, DTMs are often used from photogrammetry or satellite interferometry, such as the SRTM (Shuttle Radar Topography Mission; Farr et al., 2007) or ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer; Yamaguchi et al., 1998), freely available at https://earthexplorer.usgs.gov/ (last access: 16 October 2020). The SRTM has the widest cover and is the most commonly used terrain data set. Its latest version (V3.0 1”, 2014) is 1 arcsecond on most of the planet's surface, i.e. about 27 m×31 m resolution and an absolute height error equal to 6.2 m at Perdigão's latitude (Farr et al., 2007). With the advent of high-resolution techniques such as lidar aerial survey, terrain data have become available with resolutions above 10 m and vertical accuracy typically below 0.2 m (Hawker et al., 2018), and the question is whether such high resolution is needed in the computational modelling of atmospheric flows over complex terrain.

1.1 Literature review

Grid-independent calculations are a concept very dear to computational fluid dynamics practitioners (e.g. Roache, 1998) as a means for reducing discretization errors. However, its application in the context of atmospheric flows is not that simple because every level of grid refinement brings another level of surface detail; see for instance the coastline paradox (Mandelbrot, 1967, 1982). In this case, because the flow is driven by topography, the flow model results are directly correlated to the terrain data, and our problem is common to what can be encountered in geomorphology, with applications in hydrology (e.g. Zhang and Montgomery, 1994; Wise, 2000; Deng et al., 2007; Savage et al., 2016), where the DTM grid size affects the drainage area. In spite of its importance, to our knowledge, there is no systematic study on the appropriate grid size for resolving the terrain in microscale modelling of atmospheric flow over complex terrain.

Work has been done on quantifying the impact of using different DTMs and resolution on terrain attributes, such as elevation, slope, plan and profile curvature, and topographic wetness index. For instance, Mahalingam and Olsen (2016) note that DEMs are often obtained and resampled without considering the influence of their source and data collection method. Finer meshes do not necessarily mean higher accuracy in prediction (with examples for landslide mapping where terrain slope has a great influence), with the DEM source being an important consideration.

DeWitt et al. (2015) compared several DEMs (USGS, SRTM, a statewide photogrammetric DEM and ASTER) to a high-accuracy lidar DEM to assess their differences in rugged topography through elevation, basic descriptive statistics and histograms. Root mean square error ranged from 3 (using a photogrammetric DEM) to approximately 15 (using the SRTM) or 17 m (using ASTER).

Deng et al. (2007) indicated that the mesh resolution can change not only terrain attributes in specific points but also the topographic meaning of attributes at each point. They concluded that variation in terrain attributes was consistent with resolution change and that the response patterns were dependent on the landform classes of the area. Deng et al. (2007) introduced the concept of threshold resolution, i.e. the resolution beyond which the model quality deteriorated quickly but below which no significant improvement in modelling results was observed.

Florinsky and Kuryakova (2000) developed an experimental three-step statistical method to determine an adequate DEM resolution to represent topographic variables and landscape properties at a microscale (exemplified by soil moisture) by performing a set of correlation analyses between resolutions.

Diebold et al. (2013) showed the effect of grid size in large-eddy simulation (LES) of flow over Bolund. Lange et al. (2017) addressed the question of how to represent the small topographic features of Bolund in wind tunnel modelling, comparing a round and a sharp edge of a cliff in a wind tunnel and concluding that the cliff with the sharp edge gives an annual energy production of a wind turbine near the escarpment that is 20 % to 51 % of the round-edge case.

1.2 Objectives and outline

The objective of the present study is to determine the terrain resolution required to accurately resolve the atmospheric flow over Perdigão and mountainous terrain in general. One needs to assess the terrain horizontal resolution before assessing the effect of other (also important) causes of differences between experimental and computational results. Many computational studies based on Perdigão data are expected, and it is important to assess the terrain resolution requirements first.

In what follows, we describe the techniques used for aerial and terrestrial surveying (Sect. 2) plus the post-processing of those data and the determination of the main geometrical parameters of the Perdigão site (Sect. 3). The results of terrain attributes and computational flow modelling are the subject of Sects. 4 and 5. The paper ends (Sect. 6) with conclusions and recommendations on the most appropriate DTM and grid size required in the computational modelling of the flow over Perdigão.

2 Topographical surveying: equipment and techniques

2.1 Airborne laser scanning (2015)

The lidar aerial survey (Mallet and Bretar, 2009) was performed on 15 March 2015 by NIRAS (2015), with assistance from Blom TopEye. The survey covered an area equal to 22 071 075 m² (Fig. 1) and was completed in one session at an altitude of 500 m with a TopEye system S/N 444 and a camera mounted on a helicopter. The number of points of the lidar point cloud was approximately 993 198 375, an average point density inside the project area equal to 45 and 12.6 points m⁻² if restricted to the ground-class points (Fig. 1a and Sect. 3.1). The photography (a total of 744 photos, stored as 300 m ×300 m tiles) was performed with a Phase One iXA180 medium-format camera with 10328 pixel ×7760 pixel sensor resolution, yielding a ground resolution equal to 4.7 cm (Fig. 1b).

Lidar data were checked for point density control by Blom's software TPDS (TopEye Point Density and Statistics), the area being fully covered by lidar data with exceptions for watersheds. GPS signal was processed using data from three Portuguese reference network stations (CBRA, MELR and PORT; cf. DGT, 2017) after assistance by the Portuguese National Mapping Agency (Direção-Geral do Território, Divisão de Geodesia). Discrepancies between flight lines (based on Blom’s software TASQ, TopEye Area Statistics and Qualities, calculated on sub-areas of 1 m and after matching of 204 104 275 observations) showed a maximum altitude deviation and RMSE equal to 0.490 and 0.061 m. In 75 % of the sub-areas, the RMSE was lower than 60 mm.

Figure 1 Lidar aerial survey and orthophoto: (a) ground-class point distribution; (b) orthophoto (houses in red).

The raw data of the NIRAS (2015) campaign comprised the lidar point cloud in LASer (LAS; version 1.2) format and the orthophotos in 20 and 5 cm resolution; for more information and availability on these data see Palma et al. (2020e). The production of the digital terrain model based on the lidar point cloud is the subject of Sect. 3.1.

2.2 Terrestrial surveying (2017 and 2018)

During the installation of scientific equipment (Nov 2016–May 2017), terrain elevation was measured in situ (Palma et al., 2018) for an accurate and final determination of the elevation data of part of the instrumentation. The measuring equipment was a Leica system comprised of the following units: (1) Leica Nova Multistation MS50 (http://w3.leica-geosystems.com/downloads123/zz/tps/nova_ms50/white-tech-paper/Leica_Nova_MS50_TPA_en.pdf, last access: 20 October 2020), (2) http://leica-geosystems.com/products/gnss-systems/smart-antennas/leica-viva-gs14 (Leica Viva GS14 – GNSS Smart Antenna), (3) http://leicatotalstation.org/tag/ctr16/ (Leica CRT16 Bluetooth Cap, last access: 21 October 2020) and (4) http://leicatotalstation.org/tag/grz121/ (Leica GRZ121 360 Reflector PRO Surveying Prism, last access: 20 October 2020).

In 2017, a piece of land required changes for installation of tower 20/tse04. The topographic survey of that region was carried out (Alves, 2018) and incorporated in the lidar-based terrain model of March 2015. This survey was performed by Spectra Physics (SP60 GNNS receiver and data collector T41 with Survey Pro software) equipment and software by Sierrasoft (PROST Premium/Topko Standart, Version 14.3).

Figure 2 Workflow diagram for producing the terrain map using LAStools©.

3 Terrain model

3.1 Lidar point cloud processing

The lidar data was classified in four data type classes (ground, vegetation, unassigned and noise), stored in LAS file format, and then post-processed with tools pertaining to the LAStools^© software suite (LAStools, 2019) in three stages (Fig. 2).

Stage 1 was concerned with the extraction of the lidar raw data. The ground-class point cloud had irregular spacing (Fig. 1a), with lower point density in regions of vegetation clumping or non-overlapping scans. Some small, distinctive areas were found to be devoid of ground points due to watersheds and lidar reading or classification errors.

Stage 2 involved the reclassification of abnormal data. A first procedure was used to reclassify a particular area of noise-classified points into the ground or vegetation classes, which would otherwise be void of ground points. Points with excessive (>700 m) or negative elevations above sea level (a.s.l.) were removed during this stage. Isolated points above or below the more spatially dense point cloud, classified as ground or vegetation, were identified and removed using the lasnoise tool (LAStools, 2019).

In Stage 3, a triangulated irregular network (TIN), based on the Delaunay triangulation, was obtained for the ground-classified points. The DTM was obtained by interpolating the heights into a regular mesh with a resolution of 2 m ×2 m, the highest horizontal resolution of the terrain elevation within the Perdigão site.

3.1.1 Buildings

It is not clear whether a DTM should comprise buildings and other human-made artefacts that are usually part of digital surface models (DSMs). In the context of the present study, buildings are long-standing structures as a terrain feature, and we saw no reason for buildings to not be part of the DTM. The houses, in Fig. 1b, are family houses of about 15 m ×15 m and 5 m height that will show on the finest mesh only and as a point elevation. Unless there are a few neighbouring buildings, the ability to resolve isolated houses is limited.

The first task, to include the building data in the DTM, embraced the digitalization of all buildings from the orthophotos. This process was needed to retain the building polygons as close as possible to their exact shape and location. The second task involved the extraction of unassigned data points, which included buildings and adjacent vegetation among other structures that fell within the polygons. These points were reclassified using the lasclassify tool (LAStools, 2019) to further remove adjacent and overhanging vegetation, and the resulting building-class points were converted to the ground class. The third and final stage comprised the generation of a TIN from the new ground point cloud followed by interpolation (blast2dem tool; LAStools, 2019) of the heights to a regular mesh with a resolution of 2 m ×2 m.

Calculations including the building data showed minor or no visible flow changes and were discarded. Nevertheless, for future use, two DTM versions, with and without buildings, are made available (Palma et al., 2020e).

3.2 Two-dimensionality and main geometrical parameters

One of the reasons why Perdigão was selected was its geometry, namely the parallelism between the two ridges and their large length relative to the width, bringing the orography close to two-dimensionality.

3.2.1 Area of interest (AOI), reference lines and locations

For scaling and dimensional analysis, the main geometrical parameters of the Perdigão site were determined, and the area of interest (AOI) was defined (Fig. 3 and Table 1): a rectangular shape of approximately 3 km ×4 km, with the lower left corner at 4394131N, 608589E and aligned with the centreline (ℓ_C; SE–NW direction, 135^∘). This area, centred near station 131, included the SW and the NE ridge, the valley, and the location of most of the instrumentation deployed in Perdigão. Note that the coordinate system was converted from ETRS89 PT-TM06 (original source) to ED50 UTM29 and is used throughout the document as eastings and northings. Station number (No.) and code are given as on the Perdigão web site (https://perdigao.fe.up.pt, last access: 20 October 2020).

Figure 3 Area of interest and transects (ED50 UTM 29N).

Table 1 Reference points as in Fig. 3 (ED50 UTM 29N).

3.2.2 Terrain profile and slope

The terrain profile (Fig. 4) is not uniform along the valley, which becomes narrower and deeper along the SE–NW direction. For instance, the NW ridge height relative to a reference height (h_ref; the mean height of the surrounding terrain in a 20 km ×20 km area) equal to 250 m a.s.l. varies between 201.1 and 251.4 m, and the distance between ridges is 1358.0 and 1480.0 m on transects A and D (Table 2).

Apart from ℓ_C, six additional lines were defined: ℓ_SW and ℓ_NE along the SW and the NE ridges and A, B, C and D perpendicular to the ridges (SW–NE direction, 225^∘) and related to four main transects – A and D, which delimit the northernmost (station 39) and southernmost (station 32) locations of the great majority of the instrumentation, and transects B and C, which delimit a narrower region, determined by locations of stations 105/LRWS#5 and 20/tse04.

Other geometric variables (Fig. 4) are the height of ridges (h_SW and h_NE) and the valley (h_val) relative to the reference height (h_ref), the half-widths of the ridges ($l_{\mathrm{SWw}},l_{\mathrm{SWe}},l_{\mathrm{NEw}}$ and l_NEe) at half-height (h_SW∕2 and h_NE∕2), and the distance between ridges ℓ.

Figure 4 Terrain profile: geometrical parameters.

Table 2 Main geometric variables of transects A, B, C and D (slope S in degrees, ^∘, and length and elevation ℓ and h in metres; h_ref=250 m a.s.l.).

The ridges' orientation was determined by a linear regression of two z maxima for each j (mesh oriented with a SW–NE direction) on a 20 m grid between transects A and D (≈1650 m) and between transects B and C (≈530 m). The deviations from parallelism are 4.3^∘ if restricted to the region between A and D. Between transects B and C, where the core of the instrumentation was, the ridges were parallel within 2.8^∘, i.e. 139.1 and 136.3^∘ in the case of SW and NE ridges. The slope ($S=\left|\mathrm{atan}\right(h_{\mathrm{SW},\mathrm{NE}}/\mathrm{2}\left)\right|/{\mathrm{\ell }}_{\mathrm{SW},\mathrm{NE}}$), also on a 20 m grid, varies between 21.08 and 45.09^∘, always above the threshold for flow separation under neutral conditions (Wood, 1995).

4 Digital terrain model: results and discussion

The terrain elevation and slope are the two main terrain attributes for classification of terrain complexity and the ones to replicate accurately by terrain models. In this section, three DTMs of the Perdigão site are analysed within the AOI by comparing terrain elevation and slope on meshes based on these terrain models, with the terrain elevation and slope measured by the lidar aerial survey data within the AOI.

The three DTMs (Fig. 5 and Sect. 4.4) were the following: (1) ALS, the area sampled by lidar with a 2 m resolution; (2) military (Mil), the Portuguese Army cartography around Perdigão, 10 m horizontal resolution (available from the Portuguese Army Geospatial Information Centre; CIGeoE Centro de Informação Geoespacial do Exército, https://www.igeoe.pt (last access: 20 October 2020); a total of eight sheets – numbers 290.4, 291.3, 302.2, 302.4, 303.1, 303.3, 313.2 and 314.1 – at a scale equal to 1:25 000); and (3) SRTM, with a resolution of about 27 m×31 m. Information on availability of these data can be found in Palma et al. (2020e) and in the data availability section at the end of the present study.

Because the ALS was the highest-resolution map and the most accurate representation of the terrain in Perdigão, it was the one against which the accuracy of alternative terrain data sources was evaluated.

Figure 5 Total area comprised of SRTM, military and airborne data.

Concerning the terrestrial surveys in 2017 and 2018 (Sect. 2.2), a sample of those measurements confirmed the high quality of lidar airborne measurements. The survey carried out in 2018 showed that the terrain change due to installation of tower 20 yielded alteration of the terrain that was not significant (less than 1 m).

4.1 Mesh generation

For comparison of terrain attributes, elevation and slope, regularly spaced meshes of 80, 40, 20 and 10 m (sizes n_i×nj of 39×51, 77×101, 153×201 and 305×401, respectively) were generated within the AOI. The resampling procedure was similar to Deng et al. (2007); i.e. one out of two points was retained to assure that every point in the coarser resolutions existed in the finer ones. Coarser meshes are resampled versions of the 2 m resolution mesh, obtained by removing additional nodes.

4.2 Elevation at five reference points

Five points (Table 1) were selected for DTM comparison: towers 20/tse04, 25/tse09 and 29/tse13 – the three 100 m meteorological towers, comprising a transect aligned with the dominant wind direction – tower 37/rsw06, and the wind turbine location along the SW ridge.

Figure 6 shows the absolute error ($z_{\mathrm{error}}=z_{\mathrm{80},\mathrm{40},\mathrm{20},\mathrm{10}}-z_{\mathrm{2}}$), difference between the elevation at a given mesh and DTM source with respect to the terrain elevation on the reference mesh (ALS₂). In the case of SRTM-based meshes (Fig. 6a), the error tends to a plateau at resolutions equal to 40 m. Similar behaviour is found in the case of the Mil database (Fig. 6b) but at 20 m resolution; 20 and 10 m meshes increase the error at 20/tse04, and meshes at higher resolution than the uncertainty of this database must be avoided. Contrary to the SRTM and Mil, when using ALS (Fig. 6c) with mesh refinement there is a noticeable error reduction at all five points.

Figure 6 Impact of mesh resolution on reference points: (a) SRTM, (b) Mil, (c) ALS.

Figure 7 Impact of mesh resolution on RMSE.

4.3 Elevation and slope in the area of interest (AOI)

As the DTM quality increased from the SRTM to Mil and ALS, the maximum terrain elevation (z_Max) also increased, from 530.2 to 540.5 and 540.8 m, and the minimum (z_Min) decreased from 165.0 to 158.8 and 156.8 m (Table 3). Maxima and minima terrain elevations are set by the DTM source; maxima and minima are similar for a given DTM regardless of the grid size, which was a consequence of the procedure for mesh refinement. The 10 m difference between SRTM and the ALS values is consistent with the RMSE of the SRTM, equal to 6.2 m for Eurasia (Table 1, Farr et al., 2007), and also with the conclusions by DeWitt et al. (2015).

Table 3 Maxima and minima terrain elevation, based on SRTM, Mil and ALS data.

Figure 8 Elevation error of resampled meshes (ALS terrain data) over the AOI surface: (a) 80 m mesh, (b) 40 m mesh.

Table 4 Maxima and minima slope in the x (SW–NE, 225^∘) direction, based on SRTM, Mil and ALS terrain data.

The error distribution (Figs. 8 and 10a) shows an overprediction of the terrain elevation along the valley and an underprediction along the ridges, with both much reduced between the 80 and the 40 m resolution meshes, with the latter showing a mostly uniform error distribution of around 1 m (Fig. 8b).

The RMSE error (Fig. 7) over the whole AOI is consistent with the elevation error and the inherent uncertainty of every DTM source; with mesh refinement every DTM tends to its resolution level. The minimum RMSEs of the SRTM, Mil and ALS are 7.43, 4.66 and 0.61 m at resolutions of 10 m.

Figure 9 Slope in x (SW–NE, 225^∘) direction with different resolutions mapped on the AOI's surface (ALS terrain data): (a) 80 m mesh, (b) 10 m mesh.

Figure 10 Terrain along transect C (see Fig. 3), ALS-based meshes: (a) terrain profile and elevation error, (b) slope in the x direction.

The maximum slope (55.85 and 64.76^∘) was about 50 % higher on 20 and 10 m meshes compared with the coarser resolution (37.31 and 44.24^∘ meshes, 80 and 40 m; Table 4). The negative slope increased from −37.33^∘ to −67.81^∘ as the resolution increased from 80 to 10 m. The histogram of the slope in the x (SW–NE, 225^∘) direction (Silva, 2018) shifted to the right as the content at low slopes decreases, and more and more higher slope locations were resolved. Because the ridges are quasi-two-dimensional, the y (SE–NW, 135^∘) direction slope was residual (Silva, 2018) compared to the x direction slope (Fig. 9) and is not shown here.

The larger errors occurred at locations of higher slope (Fig. 10), and these are the locations where the grid refinement was also the most effective in reducing the elevation error. For instance, errors equal to 11.5 and −15.8 m (at $x=-\mathrm{720}$ m and x=766 m) on an 80 m mesh were reduced to 7.5 and −2.5 m on a 40 m mesh.

4.4 Spectra analysis

Spectra of terrain elevation show the ALS resolution 1 order of magnitude higher compared to SRTM data (Fig. 11). The figure also displays two scaling ranges, typical of global topographies (e.g. Nikora and Goring, 2004), with exponents equal to $-\mathrm{7}/\mathrm{4}$ and $-\mathrm{11}/\mathrm{3}$.

Figure 11 Spectra analysis for DTM and meshes: (a) SRTM, (b)  Mil, (c) ALS.

As expected, there is an increase in resolved spectral range with mesh refinement and an overlap between meshes with ALS data. In the case of SRTM- and Mil-based meshes (Figs. 11a and b), linear refinements (20 and 10 m meshes) cannot replicate the decay for higher frequencies and overcome the inherent resolution of the original data. Mesh quality was as good as the terrain data source. Only meshes based on the ALS (Fig. 11c) have the ability to reproduce accurately the high-frequency range (10⁻¹ rad m${}^{-\mathrm{1}}<k<\mathrm{1}$ rad m⁻¹). The SRTM is restricted to around 30 m resolution, and meshes 20 m ×20 m and 10 m ×10 m, with identical z_Max and z_Min (531 and 165 m), are unable to replicate the ALS measured values, z_Max and z_Min (Table 3). Grid refinement cannot go beyond the inherent resolution of the DTM.

4.5 RIX (ruggedness index)

The RIX (ruggedness index) is one of the major parameters in WAsP (Mortensen et al., 2004). It has the goal of quantifying the terrain complexity. The operational envelope of WAsP corresponds to a RIX value of 0 %. The ALS data shows a maximum of 23.7 % and an overall higher value of RIX (average 15.22 %), while the SRTM reaches 19.7 % (average 11.6 %). Lower-resolution terrain data underestimate the terrain complexity.

5 Flow modelling

In this section, the results of computational runs on different meshes are discussed. A set of experimental data (UCAR/NCAR-EOL, 2019; 30 min averaged) on 4 May 2017, 22:09–22:39 UTM, is also included. This was the day and the time when the assumption of conditions of stationarity based on measurements at tower 20/tse04 were valid according to Carvalho (2019), and the flow was non-stratified based on a bulk Richardson number (R_B) equal to −0.03

$$\begin{array}{}\text{(1)}& R_{\mathrm{B}}={\displaystyle \frac{g\left({\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{100}}-{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{2}}\right)\mathrm{\Delta }z}{{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{100}}\left[(U_{\mathrm{100}}{)}^{\mathrm{2}}+(V_{\mathrm{100}}{)}^{\mathrm{2}}\right]}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{100}}$ and ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{2}}$ are the mean potential temperature at 100 m and 2 m a.g.l. (above the ground level), Δz=100 m, and U₁₀₀ and V₁₀₀ are the mean horizontal components of the velocity vector also at 100 m a.g.l. The temperature obtained from measurement data was converted into potential temperature using the following approximation (Stull, 1988):

$$\begin{array}{}\text{(2)}& \stackrel{\mathrm{‾}}{\mathit{\theta }}\approx \stackrel{\mathrm{‾}}{T}+\left({\displaystyle \frac{g}{c_p}}\right)z.\end{array}$$

The data set choice was conditioned by the computational flow model being used. Because computational results do not consider, for instance, surface cover heterogeneity, discrepancies are expected when compared with experimental data, which are included here for guidance only.

5.1 Computational flow model

The computational code VENTOS^®/2 (cf., Castro et al., 2003; Palma et al., 2008), developed for atmospheric flows over complex terrain, was used in steady-state formulation. It solves the Reynolds-averaged Navier–Stokes set of equations for a turbulent flow (k−ε model), with a terrain-following structured mesh, also allowing the simulation of forested terrain (Costa et al., 2006) and wind turbine wakes (Gomes et al., 2014; Gomes and Palma, 2016).

5.1.1 Integration domain and boundary conditions

The model topography (domain size: 19 km ×18.8 km around the central location 4396621N, 608250E) was obtained by bilinear interpolation of terrain data. The positioning of the domain boundaries and their impact on flow variables were part of the work of Silva (2018).

At the inlet a log-law profile was set. To ensure an equilibrium shear stress, the k profile decreases with the square of height above ground level. At the top of the domain a zero shear stress condition was used. The inlet profile's development is capped at the boundary layer's limit, all quantities being constant above that height. At the lateral boundaries a symmetry condition was applied.

The ground was modelled as a rough surface, a wall function, a log-law defining the velocity at the node closest to the ground, and the turbulence model quantities k and ε. The values used in the computational model for z₀ (roughness length) and u_* (friction velocity) were 0.1 (indicated by Wagner et al., 2019, as the roughness length near the double ridge area after conversion from the CORINE Land Cover classes) and 0.25. These values were uniform for the whole domain. The surface cover (forest patches and height) and its representation in the computational model (roughness length, leaf area index, use of a canopy model) are still a work in progress as the presence of eucalyptus and pine tree patches in the area are expected to have an impact on the flow. However, this would increase the number of variables influencing the flow results, masking the effects of the digital terrain model alone. See for instance the effects of forest resolution and wind orientation relative to the forest stands in the computational modelling of flow over forests in Costa et al. (2006).

5.2 Computational meshes

A total of 18 computational meshes (Table 5) were used. The central part of the domain (4 km×6 km, based on ALS and Mil terrain data) was resolved with uniform horizontal resolution (20 m ×20 m, 40 m ×40 m and 80 m ×80 m), expanding towards the domain boundaries with factors f_x and f_y close to 1, to minimize the discretization errors. The domain's height (3000 m) was discretized by 100 nodes (N_k), with the first node 2 m above ground level and a grid expansion f_z=1.0435, yielding a maximum cell size (Δ_z) equal to 124 m. For availability of meshes ALS.SW.## and ALS.NE.##, see Palma et al. (2020e) and information in the data availability section at the end of this study.

A preliminary analysis showed that the flow variables had low sensitivity to the number of nodes in the vertical (n_k) opposed to the height of the first node above ground level, which showed a significant impact, and is worthy of further studies.

Table 5 Computational meshes (Δz_Min=2 m).

^* Did not meet residual criteria.

Three types of meshes were used: the SRTM, with the whole domain based on the SRTM data; Mil, a combination of SRTM and military maps; and ALS, based on all three DTM sources (Fig. 5). A minimum of 8818 iterations and 3.87 h of computing time and, in the case of mesh Mil.NE.20, a maximum of 20 033 iterations and 50 times more computing time were required. Number of iterations is a better indicator of the actual computing time since the value given here was influenced by the computer load at the time of the calculations.

5.3 Flow pattern

The flow modelling analysis was based on the flow patterns at two transects in the case of SW and NE winds (Figs. 12 and 13) and wind speed, wind direction and turbulent kinetic energy results for SW (Figs. 14–16).

As expected, the flow pattern (Figs. 12 and 13) is characterized by separated flow regions on the lee side of either ridge. The figures are coloured by the spanwise velocity component (v), showing two different streams: up-valley on the lee side of the SW ridge and down-valley on the upwind side of the NE ridge (Fig. 12) and down-valley in the case of NE winds (Fig. 13). The ridge height increases with the grid resolution (see insets), and the detachment point moves to higher elevations, yielding longer and deeper separated flow regions (see for instance Fig. 13).

Figure 12 Impact of mesh resolution on separation zone in transect that crosses tower 20/tse04 in the case of SW winds.

Figure 13 Impact of mesh resolution on separation zone in transect that crosses tower 29/tse13 in the case of NE winds.

Menke et al. (2019), in their analysis of the experimental data, reported average length and depth equal to 697 and 157 m for both SW and NE wind directions and stratification levels based on the gradient Richardson number between −1 and 1; in the case of neutral flow, length and height equal to 807 and 192 m were reported for a 10 min period. The length and height of the separation zone, in Table 6, tend to increase with the grid refinement (with the exception of the SW winds when refining from 40 to 20 m resolution), predicting a recirculation region longer and narrower compared with the measurements.

Table 6 Length and maximum depth of separation zone.

5.4 Southwesterly winds

The wind speed, wind direction and turbulent kinetic energy profiles at towers 20/tse04, 25/tse09 and 29/tse13 (Figs. 14, 15 and 16) show a good agreement of the wind direction with the measurements, a poor agreement of the wind speed (underprediction) at all towers, and underprediction of the turbulent kinetic energy in the valley and at the NE ridge (towers 25/tse09 and 29/tse13; Figs. 15c and 16c). A good agreement between computational and experimental results is not expected, mainly because of the uniform roughness length; the important issue is the sensitivity of the computational results to the different numerical meshes.

As an indicator of stationarity, the mean values of the experimental results over the 30 min period are plotted, showing the minima and maxima within that period as error bars. Departure from stationarity condition reaches a higher magnitude in the valley (Fig. 15), given the location of tower 25/tse09 inside the recirculation zone. Unsteadiness is a well-known characteristic of recirculation zones, and their prediction is very sensitive to spatial resolution (Castro et al., 2003) and terrain model as shown by Fig. 15b. The separated flow region, tower 25/tse09 (Fig. 15), is characterized by low wind speed and rotation of the wind with the distance above the ground. The wind speeds at z_asl>100 decrease as the mesh is refined, and the height of the recirculation zone increases (Table 6). The flow in the valley is aligned with the valley and therefore perpendicular to the ridges and the incoming wind. The predicted wind direction is in close agreement with the measurements, with the exception of 40 and 20 m meshes based on the Mil DTM.

As a whole, the flow results appear to be more sensitive to the resolution than to the DTM (see, for instance, Fig. 14a), with the results on 40 and 20 m meshes detached from results on the 80 m mesh. A resolution of at least 40 m is required.

Figure 14 Simulation results and experimental data in tower 20/tse04 for SW winds: (a) wind speed, (b) wind direction, (c) turbulent kinetic energy.

Figure 15 Simulation results and experimental data in tower 25/tse09 for SW winds: (a) wind speed, (b) wind direction, (c) turbulent kinetic energy.

Figure 16 Simulation results and experimental data in tower 29/tse13 for SW winds: (a) wind speed, (b) wind direction, (c) turbulent kinetic energy.

The profiles (not shown) in the case of NE winds (45^∘) are similar to SW winds, apart from the situation being reversed, since in this case the first and second ridge are the NE and SW ridges.

Differences between the profiles and the reference profile ALS₂₀ were measured in terms of RMSE (Tables 7, 8 and 9), which, in general, show a pattern similar to the slope (Table 4), where the RMSE decreases by either refining the mesh or, for a given mesh, moving from the SRTM to Mil- and ALS-based meshes. RMSE values at towers 20/tse04 and 29/tse13 (Tables 7 and 9) on the hills depend on the dominant wind directions, showing the effects of the valley flow on the downstream hill. The effect of calculations on 20 m mesh compared to those on 40 m mesh are less than the effect of calculations on 40 m mesh compared to those on 80 m mesh.

Table 7 RMSE of wind speed, wind direction and turbulent kinetic energy for tower 20/tse04.

Table 8 RMSE of wind speed, wind direction and turbulent kinetic energy for tower 25/tse09.

Table 9 RMSE of wind speed, wind direction and turbulent kinetic energy for tower 29/tse13.

6 Discussion and conclusions

Meshes for computational modelling of flow over the Perdigão site were created based on three digital terrain models: high-resolution (2 m resolution) airborne lidar survey (ALS), military (10 m) and SRTM (30 m) data. The mesh appraisal was carried out in two ways: by their ability to replicate the two main terrain attributes, elevation and slope, and by their effect on the wind flow computational results (wind speed, wind direction and turbulence kinetic energy) at three locations.

Regarding the digital terrain models, the main conclusions were the following:

1. The SRTM data are not an accurate representation of the Perdigão site.

2. Only meshes based on the ALS have the ability to reproduce the smaller scales between 10 and 100 m.

3. The ALS data yielded the lowest elevation errors: average RMSE around 5.8 m on 80 m, decreasing to 0.6 m on 10 m mesh.

4. The RMSE for the SRTM does not go below 7.4 m. A 40 m horizontal resolution based on the ALS data is enough to achieve an error below 1.4 m in five key locations and below 0.28 m using a 20 m mesh.

5. The maximum terrain slope was about 1.8 times higher (−67.81^∘) on a 20 m mesh resolution compared with an 80 m mesh resolution (−37.33^∘). An 80 m mesh does not accurately represent elevation and slope, mainly near the extreme elevation values (highs and lows).

The effect of the terrain model on the wind speed, wind direction and turbulent kinetic energy at three locations (SW ridge, valley and NW ridge) and two incoming wind directions (SW and NE) was the following:

1. In the case of SW winds, the mesh resolution effects on the SW ridge were restricted to the first 100 m a.g.l., where mesh refinement decreased the wind speed and degraded the quantitative agreement with the experimental data, though replicating the profile shape.

2. Separated flow field in the valley is perpendicular to the main flow direction. This region increases in height and length with the mesh refinement.

3. The flow (mainly the wind direction) in the valley was the most affected by terrain resolution; low velocities (about 1 m s⁻¹) are associated with large variations in wind direction within the first 150 m a.g.l.

Concerning the digital terrain models and meshes, the conclusions were the following:

1. It was found that 40 and 20 m meshes are resolutions – threshold resolution – beyond which no or insignificant changes occur in terrain attributes, elevation and slope, and in the flow field variables, wind speed, wind direction and turbulent kinetic energy.

2. It is recommended that at least 40 and 20 m meshes based on military and ALS be used to describe the Perdigão site, with the SRTM restricted to far-away regions.

The conclusions hold under the conditions of the present work, namely terrain data and flow model equations and conditions. Under different conditions, further validation may be required.