Determining slope- and free air ﬂow wind proﬁles for WKB Prandtl models.

During night time, when the air close to the surface cools-down and the atmosphere becomes stable, katabatic down-slope ﬂows may occur in mountainous regions. Contrary, during day time, when the sun heats-up the air close to the surface, and the atmosphere becomes unstable, anabatic up-slope ﬂows are prevalent. Up to date, slope winds in the WKB Prandtl model are determined solely by means of specifying the height z j of the maximum turbulent jet above ground. Furthermore, depending on z j parameters K 0 and h for a height-dependent eddy thermal diﬀusivity function K H ( z ) and a parameter C determining the amplitude for the Prandtl wind speeds must be speciﬁed. They are most often estimated from height-dependent wind speeds u , including slope angle α and ﬁxed Prandtl number Pr = K M /K H with K H = K H ( z ) and K M the eddy viscosity. Having estimated those parameters, friction velocity u ∗ , friction temperature θ ∗ and sensible heat ﬂux Q H may be calculated. This article takes the reverse approach: From speciﬁed slope angle α , Pr , u ∗ , θ ∗ , Q H and speciﬁed form of the eddy thermal diﬀusivity function K H ( z ) the corresponding WKB Prandtl model is identiﬁed. Furthermore, the relationship between u ∗ and θ ∗ calculated via Monin-Obukhov similarity theory for ﬂat terrain and those for the WKB Prandtl model are investigated. Using this relationship, we give hints how our new parametrization of the WKB Prandtl model may be used to determine slope ﬂows and free air ﬂows in a micro meteorological model of an alpine valley for pollutant dispersion calculations. We illustrate our derivations by applying our algorithms for the calculation of the WKB Prandtl model to four examples with two taken from Grisogono et al. (2015).


Introduction
Dynamic pollutant dispersion simulation, (e.g., De Visscher (2014)), requires the accurate modeling of wind speeds, wind directions and sensible heat fluxes. The modeling of pollutant dispersion in alpine regions is not well understood, (e.g., Rotach and Zardi (2007)). The main problems with wind modeling in alpine valleys are the complex terrain, katabatic down-wind flows occur during the stable night boundary layer and anabatic up-hill wind flows during the unstable day time boundary layer. A further problem are winds over hills that tend to increase their velocity over the top of the hill due to increased pressure during upward flow. All of those mentioned topics have been discussed in detail in Zardi and Whitman (2012). Katabatic upward slope flow has been considered for the first time by Prandtl (1942). The original Prandtl model for slope flows assumes a steady wind field and that the mountainous slope has an infinite extension. During the last 20 years the original Prandtl model has been extended by Grisogono and Oerlemans (2000), Grisogono and Oerlemans (2001), Grisogono and Oerlemans (2002) and Grisogono et al. (2015). Those articles assume still a steady wind field but extend the Prandtl model by assuming an eddy thermal diffusivity K H (z) that is varying with orthogonal height z above ground, (e.g., Grisogono and Oerlemans (2000)), resulting in the WKB-Prandtl model. Later, Grisogono et al. (2015) extend this model by means of adding to the steady state Boussinesq, (Boussinesq (1897)), Navier-Stokes equation a small -potential-temperature-flux-term that enhances along-slope advection of the potential temperature θ that is counteracted by parametrized turbulent mixing. We call this the WKB Prandtl model of order 1. To our knowledge, this WKB Prandtl model has been applied to the modeling of katabatic down-flow wind speeds at the Pasterze glacier, Austria, (Greuell et al. (1994), Grisogono and Oerlemans (2001)). Grisogono and Oerlemans (2001) show how the friction velocity u * , friction temperature θ * and sensible heat flux Q H can be calculated in the WKB =0 Prandtl model. Grisogono et al. (2015) but miss a derivation of u * , θ * and Q H as this was done in Grisogono and Oerlemans (2001). Our own derivations will show that for the WKB >0 Prandtl model the same formulas for u * , θ * and Q H may be used as for the WKB =0 Prandtl model but the height z j of the jet, where maximum wind speed occurs, must be taken from the WKB >0 Prandtl model. In our approach to the WKB Prandtl model it is important to have already available friction velocities u * , friction temperatures θ * and sensible heat fluxes Q H . A common procedure in boundary layer meteorology is to measure wind speeds, wind directions, temperature, humidity, air pressure, precipitation, cloud cover, illumination some meters above surface by means of meteorological stations and then to interpolate those measurements over the spatial area of investi-gation by means of techniques like inverse distance weighting, (e.g., Bhowmik and Cabral (2011)), spatial kernel methods, (e.g., Watson (1964)), kriging, (e.g., Krige (1951)), and so on. Our own interpolation of surface wind directions and speeds works separately for northing and easting wind speeds and uses inverse distance weighting. Having interpolated the wind speeds at the station heights above ground and having calculated sensible heat fluxes the investigated spatial area may be classified into different Pasquill-Gifford stability classes, (Pasquill (1961), Gifford (1961)). Thereafter the friction velocities u * and friction temperatures θ * may be calculated and the wind speeds can be interpolated vertically using the logarithmic wind profiles from similarity theory, (e.g., De Visscher (2014), Monin (1959) and Obukhov (1971)), depending uniquely on stability class, u * and θ * . Given starting values for wind speeds, u * , θ * , Q H the formulas for those parameters are iterated recursively until all those parameters stabilize. During these iterations, also northing, easting and vertical wind vectors are adjusted such that the mass conservation condition (divergence of wind equal to zero) is almost satisfied. To this, a Lagrangian optimization problem has to be solved stating that the new wind vectors deviate from the old ones as less as possible subject to having divergence almost zero. This leads to a relatively large linear equation system which has to be solved numerically by means of a successive overrelaxation method, (e.g., Varga (2002)). Sensible heat flux can be calculated from spatial latitude, sun elevation, air temperature, albedo, surface roughness length, plant cover and humidity. We calculate a first sensible heat flux using a formula given in Holtslag and van Ulden (1982), Holtslag and van Ulden (1983) and van Ulden and Holtslag (1984): where c 0 is the cloud cover, T 0 the surface temperature, r the albedo, γ is taken from Table 5, column 2 in van Ulden and Holtslag (1984) and with φ the sun elevation angle. According to De Visscher (2014) A.5.3, the albedo r is made dependent also on sun elevation and surface type (forest and grassland). During iterating through u * , θ * , wind speeds and formula (21) those parameters are successively changed until they stabilize. The question is now whether u * , θ * and Q H calculated via Monin-Obukhov similarity theory for almost flat terrain can be used also as parameters determining the Prandtl WKB model or whether they need a functional modification and how this modification must look like. If any modification is necessary we will state that.
The article is structured as follows. In section 2 the Navier-Stokes equations and their solutions for the WKB =0 and WKB >0 Prandtl model are given. Section 3 describes how parameters u * , θ * and Q H may be derived for the WBK Prandtl model . The most innovative part of this article is section 4 where it is shown how the WKB Prandtl model may be determined from knowing only α, u * , θ * , Q H , Prandtl number P r, and parametrically specified height dependent eddy thermal diffusivity K H (z). Section 5 makes use of some illustrating examples and shows how the WKB Prandtl model may be estimated accurately by means of knowing only the slope angle α, the friction parameters and the sensible heat flux. Furthermore, hints are given how this new parametrization of the WKB Prandtl model can be used for modeling the micro-meteorology in pollutant dispersion simulation in an alpine valley. The question of the relationship of the friction parameters u * and θ * in flat and steep terrain will be answered in section 6. In section 7 we combine the WKB Prandtl model for steep terrain with the one for flat terrain such that the complete modeling of micro meteorology in e.g. an alpine valley becomes possible. The final section 8 draws conclusions and gives hints for future developments.

The Prandtl-and the WKB Prandtl model
The classic Prandtl model has been invented by Prandtl (1942) and has been used by several authors, e.g., by Defant (1949), Smith (1979), Egger (1990) or Axelsen and van Dop (2009). A detailed derivation of the here to be discussed WKB Prandtl model has been given e.g., by Jurlina (2013): The governing steady-state one-dimensional equations in the tilted coordinate frame for a constant slope α, 0 < α < π 2 with valid quasi-hydrostatic approximation are given for wind speed u(z) parallel to the surface and potential temperature anomaly ∆θ(z) at height z orthogonal to the surface by Grisogono et al. (2015). Here, the parameters have the usual meaning: where g is the acceleration due to gravity, Prandtl number P r = K M /K H , eddy viscosity K M , eddy thermal diffusivity K H , absolute temperature at surface level θ 0 and ∈ (0, 1]. The boundary conditions are selected as the usual ones for the Prandtl model: ∆θ(z = z 0 ) = C, u(z = z 0 ) = 0, ∆θ(z → ∞) = 0 and u(z → ∞) = 0. Γ occurs in the assumed ordinary K-closure for the Navier-Stokes equation: where w is the wind speed orthogonal to the surface of the slope. Here, according to ordinary K-closure the potential temperature is defined as follows where z * is the vertical height above surface; The WKB approach consists of expanding wind speed u and potential temperature anomaly ∆θ into a power series of in the following way, e.g., Metivier (2003): We define in applications, like pollutant dispersion calculations, θ 0 as the absolute temperature at soil level z 0 , where z 0 is the roughness length orthogonal to surface; according to ordinary K-closure we set z * equal to z; the height coordinate z is assumed to be orthogonal to the slope; Γ is the potential temperature vertical gradient in standard atmosphere conditions and ∆θ is the anomaly of potential temperature deviation from background potential temperature. In the literature, (e.g., Garratt (1992)), it is proposed to set Γ = 0.0007Km −1 in the katabatic case and Γ = −0.0007Km −1 in the anabatic case; u is assumed to be the down-slope wind component in the katabatic case and then is assumed to be positive and is the up-slope wind component in the anabatic case and then assumed to be negative.
"The only new term in Eq. 2, as underlined in the thermodynamic equation, modulates the amount of nonlinearity (presumably weak in terms of regular perturbation analysis) via small parameter , where 0 ≤ 1. ", Grisogono et al. (2015). Grisogono and Oerlemans (2001) and Grisogono et al. (2015) select a heightdependent eddy thermal diffusivity function for z ≥ z 0 , as follows: which depends on three parameters related to the maximum value of K H (z) and the elevation h, where the maximum is attained, i.e., (K 0 , h, K min ). Typically, K min is given a small value, tenfold smaller than the molecular diffusivity of air.
(u 0 , ∆θ 0 ) are the solutions of the WKB =0 model and according to Grisogono and Oerlemans (2001) and Grisogono et al. (2015) are given by: where The original Prandtl solution may be recovered by means of setting K H (z) = K H . Above solution (8) retains the originality and elegance of the original Prandtl solution for K H = const and = 0: with Grisogono et al. (2015) show that first order solutions in the WKB model are given as follows: where the amplitudes are The boundary conditions for these solutions are as follows: ∆θ 1,W KB (z 0 ) = 0, u 1,W KB (z 0 ) = 0, ∆θ 1,W KB (z → ∞) = 0 and u 1,W KB (z → ∞) = 0. The total WKB solutions up to the first order then are given by: Based on empirical arguments, Grisogono et al. (2015) propose to select in the katabatic case of a stable atmosphere = 0.005 and in the anabatic case of an unstable atmosphere = 0.03. According to Metivier (2003) the first order solutions (14) and (15) are accurate up to order O( 2 ) to the two solutions (6) of the system of differential equations (2). 3 Determining u * , θ * and sensible heat flux Q H in the WKB Prandtl model Normally, the friction temperature θ * and friction velocity u * are defined as follows: where the above expressions are evaluated at height z = z 0 above ground. According to Figure 4 in Grisogono and Oerlemans (2001), the momentum flux shows a peak only some decimeters above ground and is disregarded by them. They then extrapolate the momentum flux from z j , the height where maximum wind speed occurs in the Prandtl model, linearly down to roughness length z 0 by means of calculating the tangent to the momentum flux at height z j and looking where this tangent intersects z 0 . Similarly, the temperature flux is extrapolated from z j to z 0 by assuming the value of the temperature flux below z j to be the same as that at z j , see Fig. 1. Thus, In the classical Prandtl model an explicit formula for the jet height z j may be given as follows: In the WKB ≥0 model no explicit formula for z j can be derived and must be calculated from the wind profile. According to Grisogono and Oerlemans (2001), explicit expressions for u * and θ * in the WKB =0 and the classical Prandtl model are given as follows: We have investigated numerically whether above formulas for u * and θ * from Grisogono and Oerlemans (2001) may be used also in the WKB >0 case but with height z j of the jet calculated from the WKB >0 model. To this, we have calculated numerically by means of replacing first and second derivatives occuring in formulas (18) by difference quotients, when using the WKB >0 solutions (14) and (15), and have shown that those derivatives where occurs are almost zero. The next two examples (Table 1) illustrate this fact for a katabatic and an anabatic case. Since the deviation of u * 2 from the full WKB >0 model to the one from the WKB =0 model but z j taken from the full model is very small, we propose to use Grisogono and Oerlemans (2001) formulas also for the WKB >0 model. The sensible heat flux Q H at height z j , defined as in the WKB >0 model should not be calculated from the well-known formula because the errors between them, dev. Q H in Table 1, are large. Here, i.e. c p = 1006 J/(kg · K) is the specific heat capacity of dry air at a temperature of 20 C o and ρ the density of air. Obviously, this deviation cannot be disregarded in the WKB >0 model and, as we have done, numerical approximation formulas (difference quotients) should be used to calculate Q H from formulas (14) and (21).
In using u * and θ * from formulas (20), valid for the WKB =0 model, our  Figure 3 and 6 . Calculated u * , θ * with z j from full WKB >0 model but calculated with formulas for WKB =0 model. Dev. u * 2 and dev. Q H are the errors made for non-using the full WKB >0 model for calculating u * 2 and Q H . The deviation differences are always from the full model to it's partial approximation WKB =0 .
WKB models are somewhat inconsistent with exception of the WKB =0 case.
We have tried to reproduce the sensible heat flux of the WKB model at z = z j by means of using formulas (21) and (14). But the used u * and θ * only are approximations to the true ones and actually give sensible heat flux (22). Defining we have tried to equally scale-up u * and θ * but our numerical algorithms from the next section 4 for calculating from given u * 0 , θ * 0 and Q H,0 the parameters K 0 , h and C, determining the WKB model, became numerically unstable. It seems that the relationship between u * , θ * , Q H and the parameters K 0 , h and C is very volatile and "discontinuous". Also, in our upcoming figures 3, 4, 5, 6 we have plotted only the full WKB models and have not taken into account that according to the assumptions made in Fig. 1 actually the WKB =0 model should be used for 0 ≤ z ≤ z j . According to the inconsistency in u * and θ * we could have calculated different amplitudes C for the WKB =0 part and the WKB >0 part. The C for the WKB =0 model could be defined such that the wind speeds in the WKB =0 and the WKB >0 part are identical at z j . But the potential temperature would still show a jump at z j because of the inconsistency. Thus, we decided to plot for the whole range of z only the WKB solution and to disregard the WKB =0 solution for z ≤ z j . For all that, this model still is inconsistent because (21) and (22) are different. For our application, pollutant dispersion calculations, this inconsistent model is still sufficient because, due to discretizations, we need only the average wind speeds for medium sized 3D-cubes. The only possibility to remove the mentioned inconsistencies at all would be to calculate u * and θ * numerically from (16) and (17) at jet height z = z j or roughness length z = z 0 by means of numerical derivatives. Then either the combined or the full model can be used consistently.
4 Determining the WKB Prandtl model from u * , θ * and sensible heat flux Q H .
To date, there is no explicit expression for the WKB Prandtl model known, depending only on u * , θ * and Q H . Also our approach calculates the WKB Prandtl model only numerically. We have written two R-functions, (R Core Team (2017)), for calculating the WKB Prandtl model. The first R-function calculates u * , θ * and the amplitude C of the Prandtl model when K 0 , h and Q H are given. The second R-function, Fig. 2, calls the first one in an optimization routine where the Euclidean distance between the wanted u * 0 , θ * 0 and the ones calculated from the first function is minimized. The minimization is done subject to the condition that the calculated model is proper and Q H is fixed at it's true value.
We now describe the first function: Actually, if the sensible heat flux Q H is given, equation (21) must be satisfied. If K 0 and h are given, putting to the left side of equation (21) the known sensible heat flux, we receive a quadratic equation for C. Replacing derivatives of equations (8) and (14) by difference quotients this quadratic equation can be solved analytically. It has two solutions for C and we take this solution as the appropriate one that is closest to the solution of the Prandtl WKB =0 model. Its solution also can be given analytically but without any numerical derivatives from the linear equation (21) including Q H . Starting with the WKB =0 solution for C and z j , the determination of C from the above mentioned quadratic equation is iterated. Using the jet height z j calculated (numerically) from equation (15) with the C from the previous iteration, the quadratic equation (21) for C is solved successively until convergence of C. At the end, the function calculates not only the true C but also u * and θ * from equations (20). So, what we have now from this first R-function are K 0 , h, C, u * , θ * and Q H . Grisogono and Oerlemans (2002) argue that WKB =0 models in order to be proper must fulfill the condition max(2z j , z inv ) ≤ (e 1/2 − 1)h (23) FindK0h(h l , hr, f l , fr) Return: (K 0 , h, f 0 ). where z inv is the height orthogonal to the surface where the temperature inversion takes place. Having the parameters of the WKB =0 model fixed those two parameters z j and z inv can be calculated numerically. We assume the same condition (23) to must be fulfilled also for the WKB >0 model. As already stated, we want to numerically parameterize the WKB >0 model in u * 0 , θ * 0 , Q H,0 by solving a quadratic optimization problem: with c = 0 if condition (23) is true and c = 0.1 otherwise. Here, θ * , u * are the parameters dependent on K 0 , h and Q H,0 that can be calculated by means of our first function. Thus, we have to minimize the above function in K 0 and h because according to our first function, K 0 , h and Q H,0 uniquely determine u * and θ * . Our procedure to minimize this function is in two parts. It works by numerically minimizing f in K 0 with h fixed (the L-BFGS-B -method is used, R Core Team (2017); step 2 in the function FindK0h, Fig. 2). Having calculated the optimal K 0 for several h we select the optimal tupel (K 0 , h): We use an interval bisection method. The central h in the current interval with it's suboptimal K 0 and f (K 0 , h, Q H,0 ) are compared to the corresponding f s on the left side of h and the right side of h. If the middle h and its corresponding f is the minimum of all three h (left h, middle h, right h) then the algorithm stops. If any f corresponding to h on the left or the right side of the current middle h is smaller than the middle f , then the interval is bisected and the new h is in the middle to the left end or to the right end from the middle h. The algorithm terminates when either the change in function values f is too small, < 0 , or when the f corresponding to the middle h is smaller than the left and the right f . The algorithm converges and the condition (23) is fulfilled if the function value f provided by this algorithm is smaller than 10. Then the parameterization of the WKB Prandtl model is such that the resulting u * and θ * deviate from the wanted u * 0 and θ * 0 relatively by no more than ≈ 10%. A problem with this approach are starting values K 0 for the L-BFGS-B-algorithm and a starting interval where the true h is supposed to be found. Those parameters are thus application dependent (wind field in an alpine valley, katabatic wind over a glacier,...,). We have implemented that once condition (23) is fulfilled for the middle (K 0 , h) all corresponding subintervals use this optimal K 0 -value as a starting point for the L-BFGS-B-routine. Thus, our algorithm builds-up a binary tree of suboptimal solutions (h, K 0 , f (K 0 , h, Q H,0 )) among which the one with minimal f (K 0 , h, Q H,0 ) is selected as overall solution. This binary tree is programmed as a R-function, Fig. 2, that recursively calls itself until termination.

Numerical examples
In this section we describe four examples showing how well our algorithms for calculating the new parameterization in u * , θ * and Q H work. The first two examples (table 2) are our own ones and demonstrate the parameterization for the more interesting WKB=TRUE case. The next two examples (tables 3 and 4) are taken from Figure 3 and 6 in Grisogono et al. (2015) and demonstrate both, the WKB=FALSE and the WKB=TRUE case. Each table results from three R-function calls. The first call has as input K 0 , h and C and calculates Q H as output and some other values, like u * and θ * . The second call uses K 0 , h and Q H (from the first call) to recalculate C from the first call and also determines u * and θ * . The last call uses u * , θ * and Q H (all from the first call) to recalculate K 0 , h and C. In above tables inputs are marked with (i) and outputs with (o). Outputs which may differ between the three calls are written down multiple times. Among the already mentioned calculated outputs, u * , θ * , Q H , C, K 0 and h the other outputs are u(z j ), z j , f (K 0 , h, Q H,0 ), z inv and if condition (23) is satisfied. Obviously, condition (23) is always satisfied in examples 1 and 2, but not satisfied in the WKB=TRUE case in the Grisogono examples. This is also obvious from the optimum function values f , which then are larger than 10. In all tables it is obvious from the function values f and the u * -, θ * -and Q H values that our reparameterization approach works quite well. This is also obvious from Figures 3, 4, 5 and 6 which show wind speeds and potential temperatures before and after the reparameterization corresponding to table 2, 3 and 4. Thus, in the mentioned figures only those corresponding to columns 1 and 2 in the tables give fully consistent wind speeds and potential temperatures. According to section 3 wind speeds and potential temperatures corresponding to column 3 are not fully consistent, but sufficient for our purposes. The convergence of the third R-function has some slight dependence on the starting value K 0 that is used in the L-BFGS-B optimization routine. This value should be within 50% to 150% of the true unknown K 0 -value. Otherwise, as we have observed, convergence would not occur. In the third function, there exists also a dependence on the proposed interval for the h-value which should also range from 50% to 150% of the true h-value which, as already stated, is calculated by means of bisection. The Q H -value calculated from the first function remains constant throughout all two other function calls. Normally, the call of all three functions takes less than one second. Only in certain cases where the interval for the h-value is not well specified the third function for calculating the WKB model parameterized in u * , θ * and Q H can take longer. Currently, all functions are implemented in the R-programming language with the exception of calculating the Prandtl wind speeds, z j and z inv which are implemented in C. Further time improvement of the third function would surely result when all those functions are implemented in C, too, because the third function works by calling also the second one during the interval bisection. Time improvement is important because we want to use the parameterization of the WKB Prandtl model in u * , θ * and Q H in our pollutant dispersion model where every five minutes slope winds have to be calculated from those parameters on a fine grid. In our pollutant dispersion model those parameters result from interpolation of the wind speeds of the weather stations at the height of two meters from inverse distance weighting and application of Monin-Obukhov similarity theory. Fig. 3 Potential temperature and wind speed profiles for our example 1. The output is almost identical for all three function calls. The good convergence of the algorithm can also be seen in Table 2. Especially the terminating function value f (relative error of u * and θ * ) is very small.

Fig. 4
Potential temperature and wind speed profiles for our example 2. There is a slight difference in the output between the first two (dashed lines) and the third (dotted lines) function calls. This slight difference is also visible from the terminating function value f in Table 2.
6 Translating u * and θ * for flat terrain to those for slopes and the WKB Prandtl model As already stated in the introduction we need the WKB Prandtl model for determining slope winds in the micro-meteorology of a pollutant dispersion model in an alpine valley. What we get from our model in a first step are u * , θ * , Q H determined by means of Monin-Obukhov similarity theory for flat terrain. This section will now show how those parameters can be transformed into ones for the WKB Prandtl model in steep terrain. We make the following conventions: u * , θ * are the friction velocity and friction temperature for flat terrain and u * 0 , θ * 0 are the friction velocity and friction Table 3 Input and output parameters for Fig. 3 and 6 in Grisogono, 2015 using WKB=FALSE.
Par.  Table 3 the relative error of u * and θ * is somewhat larger than in the previous examples. Table 4 Input and output parameters for Fig. 3 and Fig. 6  temperature for steep terrain with slope angle α. Furthermore, we assume that u * , θ * , Q H , u * 0 , θ * 0 and Q H,0 are defined as follows: where z 0 = z 1 cos α is the roughness length of the Prandtl model and z 1 is the roughness length for the flat terrain. Also, u, w are horizontal and vertical velocities and u 0 , w 0 are velocities along and orthogonal to the slope of a hill. Accordingly, θ and θ 0 denote vertical potential temperature and potential temperature orthogonal to the slope of the hill. Note: θ 0 (z 0 ) = θ(z 1 ).
In this section, we will prove that the relationship between above friction variables is as follows: This can be seen from the relationships described by Fig. 7 and the now following argumentation. In the proof we have to consider two cases: Downand up-slope flow.  Table 4 is larger then 10 because condition (23) is not fulfilled. For all that the relative error of u * and θ * is relatively small for Fig. 3 from Grisogono, 2015. Firstly, we consider the down-slope flow case. From the left part of Fig. 7 the following relationships are visible: β = arctan w u l = u 2 + w 2 u 0 = cos (α + β)l w 0 = sin (α + β)l.
Letting vertical height tend toward z 1 , orthogonal height above surface tends toward z 0 and, as a consequence, u and w tend both toward 0 in probability. We have to consider the limiting expressions u * 2 0 = − u 0 , w 0 z0 and θ * 0 u * 0 = − θ 0 , w 0 z0 . According to Slutsky's theorem, Slutsky (1925), we may take the limit in − u 0 , w 0 separately for u 0 and w 0 . We observe that for w converging in probability to 0, β becomes 0 and l = u. On the other hand, if u converges in probability to 0, β becomes π 2 and l = w. Putting all together we receive: sin (α + β)l = sin (α + π 2 )w = cos (α)w.
Note, that u and w are still observed at height z cos (α) vertical above ground in above formulas. Accordingly, u 0 and w 0 are observed at height z orthogonal to the slope. Thus, we have to calculate another limit to get the result of the claim: Analogously: Thus, from equations (25) and (26) it follows u * 0 = cos (α)u * , θ * 0 = θ * and Q H,0 = cos (α)Q H , which concludes the proof for the down-wind case. The proof for the up-wind case works similarly by observing the right part of Figure  7. We get: should be almost equal to the background potential temperature gradient Γ = 0.003. Our calculations show that this gradient actually is equal to 0.0512 in Fig. 8. Thus, another optimization routine would be necessary to find the optimal Γ . This could be done for example by means of adding under the square root the expression to the formulation of the optimization problem eq. (24) and optimizing additionally in Γ . Another approach would be to consider the code in Fig. 2 and the resulting minimized function value f min in equation (24) as a function of Γ and minimize with the resulting Γ equation (29). Using those approaches also the model for flat terrain u f lat , θ f lat must be redefined and recalculated such that the potential temperature θ f lat (z * ) is equal to θ 0 + Γ · z * , for z * ≤ h * 0 . These optimization problems are a topic of future research and will be considered in a forthcoming paper. Our own applications of the previously mentioned theory are concerned with the modeling of micro meteorology in an alpine valley for the purpose of pollutant dispersion calculations. The alpine valley has a length of 18km and the width from top to top of the mountains is about 6km. We have distributed 13 weather stations in the valley and at the slopes of the mountains. Those stations measure the wind speeds u 1 , ..., u 13 actually in horizontal direction. Multiplying those wind speeds by cos (α i ), where α i is the slope angle at the ith station, i = 1, ..., 13, we actually get the wind speeds parallel to the slopes at heights z i = 2 cos (α i ), where the vertical station heights are 2 meters. From those wind speeds we should actually calculate the 13 WKB Prandtl models for those stations. The results from those calculations should be 13 different u * 0 , θ * 0 and Q H,0 . We have discovered that our numerical calculations of those WKB Prandtl models just from the wind speeds are very unstable. Therefore, we have decided to follow a different approach: We assume that the wind speeds u 1 , ..., u 13 at height 2m are coming from the flat Monin-Obukhov theory logarithmic models. Using this assumption we interpolate those wind speeds at station heights 2m on a square grid separately for northing-and easting wind direction coordinates by means of inverse distance weighting. From those interpolated wind speeds we calculate by means of flat Monin-Obukhuv theory the u * , θ * and Q H which by means of equation (27) then can be transformed to the corresponding parameters for the WKB Prandtl models. Thus, our models are slightly biased. Nevertheless, they can be justified also by the uncertainty of the wind speeds themselves. According to Fig. 8, Prandtl wind speeds and Monin-Obukhov flat terrain wind speeds are almost identical at station height 2m except for a scaling factor cos (α).

Conclusions
We have demonstrated how the WKB >0 Prandtl model may be parameterized in u * , θ * , Q H and have shown several examples how well our algorithms Fig. 8 Left: WKB Prandtl model (bold line: wind speed, dotted line: potential temperature) from Figure 6, upper row, consistently combined with the standard model for flat terrain from Monin-Obukhov similarity theory. u * 0 and θ * 0 are given in Table 1. Corresponding u * and θ * can be calculated by means of equation (27). Monin-Obukhov length is calculated from L 0 = u * 2 T 0 /(0.4 · g · θ * ), where T 0 = 273.14K and g = 9.81ms −2 . The boundary layer height is determined as follows: h mix = min(0.4 u * L 0 |k| , 0.3 u * |k| ) where Lat = 46.8 • and k = 2(7.292 · 10 −5 ) sin (Lat). Those parameters uniquely determine wind speed u f lat and potential temperature θ f lat for stable atmosphere, see e.g., De Visscher (2014). Above h mix = 124.133m, both wind speed and potential temperature are assumed to be constant. Right: Relative error used in equation (28). The minimum relative error 0.02 is attained at h 0 = 15.093m, the Prandtl layer height.
work. We hope, that this new parameterization would be useful for micrometeorology and pollutant dispersion calculation where generally u * , θ * and Q H are specified from Monin-Obukhuv similarity theory for flat terrain. Problems in using our R-functions and parameterization may result when the range of K 0 -and h values is too large because then convergence of our algorithms being based on the L-BFGS-B routine and the interval bisection is not guaranteed. Another problem is condition (23) which actually should be fulfilled only for the WKB =0 model. A similar condition for the WKB >0 model is not known up to date and a topic of further research. In absence of knowing such a condition we have used condition (23) for the WKB =0 model also in the WKB >0 case. In our algorithms it is furthermore interesting that the sensible heat flux once calculated as in Holtslag and van Ulden (1982), Holtslag and van Ulden (1983), van Ulden and Holtslag (1984), remains constant and is reproduced in all our three R-functions. We have also shown how to translate u * , θ * and Q H values from Monin-Obukhov similarity theory for flat terrain to those for alpine terrain with slopes. The definition of the Prandtl layer height h 0 in equation (28) makes our theory especially useful for pollutant dispersion simulation when modeling of the complete micro meteorology (free air/slope) over complex terrain is necessary.