A Non-Linear Mixed Spectral Finite-Difference 3-D model for planetary boundary-layer flow over complex terrain

The Non-Linear Mixed Spectral Finite-Di fference (NLMSFD) model for surface boundary-layer flow over complex terrain has been extended to planetary boundary-layer flow over topography. Comparisons are made between this new version and the surface layer model. The model is also applied to simulate an Askervein experimental case. The results are discussed and compared with the observed field data.


Introduction
The Mixed Spectral Finite-Difference (MSFD) model was originally developed by Beljaars et al. (1987).It is based on the idea that the topography produces a perturbation to a steady, neutrally stratified, non-evolving flow over horizontally homogeneous flat terrain.A number of efforts have been made to improve the model calculation of the turbulent boundary-layer flow over complex terrain.Ayotte et al. (1994) evaluated the model predictions with a number of different closure schemes which range from the simple firstorder κ − Z closure to the full second-order closure.Xu and Taylor (1992) and Xu et al. (1994) made the non-linear extension of the model by including all the neglected terms.In the non-linear version of the MSFD model (NLMSFD), the model equations were solved iteratively.Another model improvement is the extension to the stable boundary layer (MSFD-STAB, see Weng et al., 1997).
Although the MSFD and NLMSFD models have been improved since the late 80's, these models can only formally apply to the surface-layer flow due to the model assumption that upwind or zero-order profiles of mean and turbulent variables are simple logarithmic surface-layer profiles, e.g., wind speed is logarithmic, shear stresses are constant and the effect of Coriolis force is absent.Ayotte and Taylor (1995) made the first effort to extend the model to the planetary boundary-layer flow with the full second-order turbulence closure (MSFD-PBL) but it was still a linear model.

The model equations
For the PBL boundary-layer flow, the model uses the Reynolds-averaged equations for steady-state, neutrally stratified incompressible flow.They are, in tensor notation including use of the summation convention, where U i and u i are the i-th component of the mean and turbulent flow respectively; f is the Coriolis parameter; i j3 is the alternating unit tensor; and angle bracket ( ) denotes an ensemble mean.The pressure gradient force consists of a nonhydrostatic mesoscale pressure component, P, and a synoptic-scale component, f i j3 U g j , where U g j is the j-th component of the geostrophic wind.To close the system of the equations, a turbulent closure scheme is needed.Weng and Taylor (2003) have shown that the so-called simple E − turbulence closure scheme performs quite well in modelling the PBL flow in most atmospheric conditions compared with more sophisticated schemes.The E − turbulence closure is a 1 1 2 -order scheme in which the prognostic equations for the turbulent kinetic energy (E) and a diagnostic equation for the turbulent length scale ( ) are used.The turbulent fluxes are locally related to mean vertical gradients and an eddy diffusivity, see Weng and Taylor (2003) for details.
Published by Copernicus Publications.

Numerical scheme and boundary conditions
As with the previous versions of MSFD/NLMSFD model, the model equations are transformed from the original standard right-hand coordinate system (x,y,z) with z in the vertical direction to a new system (X,Y,Z), by using a terrainfollowing coordinate transform.In addition, to ensure sufficient resolution near the surface and to resolve strong gradients, a log-linear coordinate transform is further used for the vertical coordinate Z.
Fourier transformation is performed in the horizontal directions to the coordinate transformed governing equations.The model variables are decomposed into an unperturbed or zero-order part, independent of x and y, corresponding to equilibrium flow over uniform flat terrain and a perturbation part due to topographic forcing.Collecting the first-order perturbation terms and solving the resulting system of equations, forms the linear version of the model (MSFD-PBL).Treating the neglected high-order terms as the source terms and solving the resulted equations iteratively leads to the NLMSFD-PBL.
After all these transformations, equations are discretized.A staggered vertical grid is used, where U-grid points are located midway between neighboring W-grid points.Variables stored at U-grid points are U, V, and P; while W and turbulent quantities are at W-grid points.The lower and upper boundaries are at W-grid locations.The resulting set of difference equations are solved using a block LU factorization algorithm (Karpik, 1988).
The surface boundary conditions used are a non-slip condition for velocity, the vertical derivative of the perturbation pressure is zero and local equilibrium condition for the turbulent quantities.At the upper boundary, the effects of the topography vanish.We set the perturbations of mean variables and the vertical derivatives of turbulent quantities to zero.The model uses periodic boundary conditions in X and Y directions.

Upstream profiles
The upstream or undisturbed profiles for the current PBL model are the results of an integration of a 1-D unsteadystate form of the model equations to quasi-steady state.For the given conditions of the site location ( f ), the geostrophic wind speed (U g , V g ) and the surface roughness length (z 0 ), we can obtain the necessary equilibrium profiles by running the 1-D PBL model of Weng and Taylor (2003).This model will be used for providing all required upstream profiles in our simulations.

Results and discussions
Model runs have been carried out for boundary-layer flow over an idealized isolated 3-D terrain and the Askervein Hill 2 WENSONG WENG AND PETER A. TAYLOR: NLMSFD FOR PLANETARY BOUNDARY-LAYER FLOW

Numerical Scheme and Boundary Conditions
As with the previous versions of MSFD/NLMSFD model, the model equations are transformed from the original standard right-hand coordinate system (x,y,z) with z in the vertical direction to a new system (X,Y,Z), by using a terrain-following coordinate transform.In addition, to ensure sufficient resolution near the surface and to resolve strong gradients, a log-linear coordinate transform is further used for the vertical coordinate Z.
Fourier transformation is performed in the horizontal directions to the coordinate transformed governing equations.The model variables are decomposed into an unperturbed or zero-order part, independent of x and y, corresponding to equilibrium flow over uniform flat terrain and a perturbation part due to topographic forcing.Collecting the first-order perturbation terms and solving the resulting system of equations, forms the linear version of the model (MSFD-PBL).Treating the neglected high-order terms as the source terms and solving the resulted equations iteratively leads to the NLMSFD-PBL.
After all these transformations, Equations are discretized.Flow variables are stored on a staggered vertical grid, where mean variables (U , V and W ) are at layer midpoints and turbulent quantities (E and turbulence fluxes) at layer lower boundary levels and z t (the top of the computation domain).The resulting set of difference equations are solved using a block LU factorization algorithm Karpik (1988).
The surface boundary conditions used are a non-slip condition for velocity, the vertical derivative of the perturbation pressure is zero and local equilibrium condition for the turbulent quantities.At the upper boundary, the effects of the topography vanish.We set the perturbations of mean variables and the vertical derivatives of turbulent quantities to zero.The model uses periodic boundary condition in X and Y directions.

Upstream Profiles
The upstream or undisturbed profiles for the current PBL model are the results of an integration of a 1-D unsteady-state form of the model equations to quasisteady state.For the given conditions of the site location (f ), the geostrophic wind speed (U g , V g ) and the surface roughness length (z 0 ), we can obtain the necessary equilibrium profiles by running the 1-D PBL model of Weng and Taylor (2003).This model will be used for providing all required upstream profiles in our simulations.

Results And Discussions
Model runs have been carried out for boundary-layer flow over an idealized isolated 3-D terrain and the  (Weng and Taylor , 2003) for the given condition of z0 = 0.03 m, f = 10 −4 s −1 , (Ug,Vg) = (13.77,−5.95)m s −1 , neutral thermal stratification.Logarithmic wind profile for the surface layer model is also inclcuded.
Askervein Hill-the site of a detailed and much referenced field study of boundary-layer flow over low hills in the 1980s.

Flow over Idealized Terrain
For the idealized case, an "isolated", "cosine-squared" terrain surface is used, which is described by where r = x 2 + y 2 and the maximum slope is πh/λ.For our test case, the values of h = 75 m and λ = 1500 m are used and the maximum slope is 0.157.Our computational domain is set as (X,Y ) ∈ [−3000,3000] m and 4000 m in the vertical and 129 × 129 × 101 grid points are employed.
Figure 1 shows the initial background profiles of U 0 , V 0 , |U 0 |, E 0 , uw 0 and vw 0 .This is the result of the 1-D PBL model runs of Weng and Taylor (2003) for surface roughness z 0 = 0.03 m, Coriolis parameter f = 10 −4 s −1 , geostrophic wind is constant with height and set to |U 0 | = 15 m s −1 (the components are selected as (U g ,V g ) = (13.77,−5.95)m s −1 so that the near surface (at z = 10 m) wind direction is approximately 0) and we assume neutral thermal stratification.This leads to the surface friction velocity, u * ≈ 0.47 m s −1 , which is used to calculate the logarithmic wind profile for the surfacelayer model runs.As can be seen clearly from the figure, the PBL has a near logarithmic mean wind profile associated with a well developed constant stress layer to Adv.Sci.Res.
www.adv-sci-res.net-the site of a detailed and much referenced field study of boundary-layer flow over low hills in the 1980s.

Flow over idealized terrain
For the idealized case, an "isolated", "cosine-squared" terrain surface is used, which is described by where r = x 2 + y 2 and the maximum slope is πh/λ.For our test case, the values of h = 75 m and λ = 1500 m are used and the maximum slope is 0.157.Our computational domain is set as (X,Y) ∈ [−3000,3000] m and 4000 m in the vertical and 129 × 129 × 101 grid points are employed.
Figure 1 shows the initial background profiles of U 0 , V 0 , |U 0 |, E 0 , uw 0 and vw 0 .This is the result of the 1-D PBL model runs of Weng and Taylor (2003) for surface roughness z 0 = 0.03 m, Coriolis parameter f = 10 −4 s −1 , geostrophic wind is constant with height and set to |U g | = 15 m s −1 (the components are selected as (U g ,V g ) = (13.77,−5.95)m s −1 so that the near surface (at z = 10 m) wind direction is approximately 0) and we assume neutral thermal stratification.This leads to the surface friction velocity, u * ≈ 0.47 m s −1 , which is used to calculate the logarithmic wind profile for the surfacelayer model runs.As can be seen clearly from the figure, the PBL has a near logarithmic mean wind profile associated with a well developed constant stress layer to a depth of about 40 m (where TKE is about 90% of its surface value).Above this surface layer, the mean wind profiles a depth of about 40 m (where TKE is about 90 % of its surface value).Above this surface layer, the mean wind profiles show a smooth blending to geostrophic values at the upper boundary layer of the model.The PBL has a depth of about 900 m.The turbulent quantities (TKE and shear stress) decrease appreciably over the depth of the boundary layer and diminish to near zero in the absence of shear at the top of the PBL.There are also supergostrophic winds around z = 650 m and the wind vector has Ekman spiral behaviour.
Figure 2 shows the comparison of vertical fractional speed-up (∆S) profiles at five different locations along the central line (y = 0) from NLMSFD-PBL and NLMSFD-SL model runs.∆S is defined as the difference between the local (U (X,Y,Z)) and the upstream wind speed (U 0 (X,Y,Z)) divided by the upstream wind speed.It can be clearly seen that the PBL version of the model predicts larger speed-up around the hill top areas and larger wind deduction at both upwind and downwind hill foot areas than the surface-layer model.Even at this lower maximum slope of 0.157, the difference between the two models is apparent.

Flow over Askervein Hill
Askervein is a 116 m high (126 m above the sea level) hill on the west coast of South Uist, one of the islands of the Outer-Hebrides (Scotland).For our model computation, the terrain data was originally prepared by Walm-sley, see Walmsley and Taylor (1996).There were a lot of very small terrain features which lead to a very noise FFT (Fast Fourier Transform) field in the original, socalled Map B. For our model runs, a 9-point smoothing was used.There is a slightly decrease in amplitude in the resulting topography.The maximum height of the hill becomes 115.73 m and the maximum slope is about 0.39, see Figure 3, where three tower lines and 10 m measurement locations are overlaid.Line-A and Line-AA go through the hill top (HT) and hill centre (CP) respectively and are perpendicular to the ridge, and Line-B goes through HT and CP along the hill ridge.
Our computational grid consists of 129 × 129 points uniformly distributed over an area of 6000 × 6000 m 2 .This makes the resolution 46.875 m and the domain sufficiently large to avoid upstream effects due to the periodic boundary conditions.The model run is for the wind direction of 210 • , -13 • from the 223 • orientation of lines A and AA.The surface roughness is assumed as uniform and taken as z 0 = 0.03 m.
To closely match the observational data, |U 0 | = 19.25 m s −1 is used.Again the 1-D PBL model of (Weng and Taylor , 2003) is used to compute the upstream profiles, which are very similar to those in Figure 1.
Model results of fractional speed-up ratio at 10 m along the line AA together with observational data are shown in Figure 4.The NLMSFD-PBL predicts slightly smaller values of ∆S over the upwind hill foot and lee www.adv-sci-res.netAdv.Sci.Res.quantities (TKE and shear stress) decrease appreciably over the depth of the boundary layer and diminish to near zero in the absence of shear at the top of the PBL.There are also supergostrophic winds around z = 650 m and the wind vector has Ekman spiral behaviour.
Figure 2 shows the comparison of vertical fractional speed-up (∆S ) profiles at five different locations along the central line (y = 0) from NLMSFD-PBL and NLMSFD-SL model runs.∆S is defined as the difference between the local (|U(X,Y,Z)|) and the upstream (|U 0 (Z)|) wind speeds divided by the upstream wind speed.It can be clearly seen that the PBL version of the model predicts larger speed-up around the hill top areas and larger wind reduction at both upwind and downwind hill foot areas than the surface-layer model.Even at this lower maximum slope of 0.157, the difference between the two models is apparent.

Flow over Askervein Hill
Askervein is a 116 m high (126 m above the sea level) hill on the west coast of South Uist, one of the islands of the Outer-Hebrides (Scotland).For our model computation, the terrain data was originally prepared by Walmsley, see Walmsley and Taylor (1996).There were a lot of very small terrain features which lead to a very noisy FFT (Fast Fourier Transform) field in the original, so-called Map B. For our model runs, a 9-point smoothing was used.There is a slightly decrease in amplitude in the resulting topography.The maximum height of the hill becomes 115.73 m and the maximum slope is about 0.39, see Fig. 3, where three tower lines and 10 m measurement locations are overlaid.Line-A and Line-AA go through the hill top (HT) and hill centre (CP) respec- tively and are perpendicular to the ridge, and Line-B goes through HT and CP along the hill ridge.Our computational grid consists of 129 × 129 points uniformly distributed over an area of 6000 × 6000 m 2 .This makes the resolution 46.875 m and the domain sufficiently large to avoid upstream effects due to the periodic boundary conditions.The model run is for the wind direction of 210     To closely match the observational data, |U g | = 19.25 m s −1 is used.Again the 1-D PBL model of (Weng and Taylor, 2003) is used to compute the upstream profiles, which are very similar to those in Fig. 1.
Model results of fractional speed-up ratio at 10 m along the line AA together with observational data are shown in Fig. 4. The NLMSFD-PBL predicts slightly smaller values of ∆S over the upwind hill foot and lee side ares and slightly larger value over hill top (HT) area and agrees better with the measured data than NLMSFD-SL results.Figure 5 shows the comparison of vertical profiles of ∆S at HT.The PBL version of the model predicts a larger value of ∆S and again gives a better agreement with the observational data than the surface-layer version of the model.

Summary
A Non-Linear Mixed Spectral Finite-Difference model for neutral planetary boundary-layer flow over complex terrain has been developed.The model uses E − turbulence closure.Some of early limitations on the surface-layer version of model are removed and model can simulate flows that are more representative of the real atmosphere.This new model has good potential in wind energy applications.

Figure 2 .
Figure 2. Comparison of vertical fractional speed-up profiles at five different locations along the central line (y = 0) of a circular cosine-squared hill from NLMSFD-PBL and NLMSFD-SL model runs.

Figure 2 .
Figure 2. Comparison of vertical fractional speed-up profiles at five different locations along the central line (y = 0) of a circular cosinesquared hill from NLMSFD-PBL and NLMSFD-SL model runs.

Figure 3 .Figure 4 .Figure 3 .
Figure 3. Computational domain of the Askervein hill and measurements twoer along the A-line, AA-line and B-line.

Figure 3 .Figure 4 .
Figure 3. Computational domain of the Askervein hill and measurements twoer along the A-line, AA-line and B-line.

Figure 5 .
Figure 5. Vertical profile of fractional speed-up ratio, ∆S, at the hill top (point HT) of Askervein hill.Comparison of model results and experimental data.

Figure 4 .Figure 5 .
Figure 4. Fractional speed-up ratio, ∆S , for flow over Askervein hill at a height of 10 m above topography along line AA.Comparison of model results and experimental data.R: NLMSFD FOR PLANETARY BOUNDARY-LAYER FLOW

Figure 5 .
Figure 5. Vertical profile of fractional speed-up ratio, ∆S , at the hill top (point HT) of Askervein hill.Comparison of model results and experimental data.
• , −13 • from the 223 • orientation of lines A and AA.The surface roughness is assumed uniform and taken as z 0 = 0.03 m.