3. Martin GERA : "Improved representation of boundary layer"
During the last period of my stay at the Royal Meteorological Institute of Belgium, supported by ALATNET Training Network, I concentrated my work on the following scientific topics:
· Investigation of a numerical scheme for solving the TKE (Turbulent Kinetic Energy) equation
· Testing some derived and suggested methods
· Comparison of chosen methods
· Results analysis
· Writing a final scientific report
One needs to solve a second order nonlinear partial differential TKE equation. It is not possible to find an analytical solution. We have to find an appropriate numerical scheme, which converges to the "true" result. Considering the vertical mesh of ALADIN one must use the finite element method for numerically solving TKE equation.
· For the linearized version of the TKE equation we found that in some cases there can exist levels where total reflection of TKE occurs.
When we analyse the energy-density function of the TKE, the reflection properties of linearized TKE depend on the wavenumbers (for long waves it can be a total barrier, but for short waves it can be a transparent level). For these reasons, the linearized version cannot be used. In some meteorological conditions the solution is in resonant mode, which abnormally increases the amplitude of the solution. However, from the linearized version of the equation we get some knowledge about the equation properties.
· In the next step we try to solve directly the TKE equation in nonlinear mode.
The purpose of solving the TKE equation, after rewriting to discrete form is to find a zero of a system of KLEV nonlinear functions in KLEV variables by a variant of the Powell hybrid method, where KLEV is the number of vertical levels in the numerical model. The main principle of this method is to find an approximate solution by iteration. The Jacobian, which is necessary for finding the gradient is calculated by a forward-difference approximation.
· Next we try to take advantage from previous methods and avoid known problems. We try to apply a Predictor-corrector method.
Here we performed a deep analysis for the choice of the linear operator. For this operator we use information form stationary TKE. Using the average values of stationary TKE as a background term, the reference state diverges from the real profile, but it is one possible choice for the linear operator. Alternatively the linearization can be done with help of a simplified TKE equation. For zero tendencies, TKE equation can be reduced to the stationary TKE equation. This equation can be linearized around the stationary TKE. Now there is no average profile of stationary TKE. But this form is more "pure" for other terms.
A next possibility is to combine the previous linear operators by averaging the dissipation term.
The solved problems are more complicated than to find stable scheme. The roots of our system should be positive, because they represent TKE or square of TKE. For this reason we expect a "balance" behaviour. Mathematically we can formulate the definition of our problem in the following way. The inverse equation operator should have positive matrix elements and the right-hand side must be positive. In this case we have positive real roots with positive right (forcing) side (this definition may be incomplete).
To have positive forcing we need to analyse the buoyant and shear terms together with the dissipation term. For the instable case there is no problem, the forcing term is always positive.
For the stable case the situation is more complicated. In the neutral case momentum and heat flux are differed by inverse Prandtl number. It means that the drag coefficient for heat flux is larger than the drag coefficient for momentum flux for neutral air stratification. The dominant contribution to the final forcing term is from the buoyant part at this moment. It means that F(z,t)<0, where F(z,t) is the forcing term which contains shear and buoyant production or consumption terms in the TKE equation.
At this moment dissipation term and source terms have an important role. The sign of F(z,t) should be changed by these terms. For these properties it is important to have restriction of the critical Richardson number. It should not have large values, so as not to have a big negative value for F(z,t). These features are significant in upper (model levels) atmosphere, where stable stratification of atmosphere is observed frequently.
In the PBL layer or in the border troposphere and tropopause other problems create obstructions. In these layers we frequently observe strong turbulent mixing. It has an influence on the mixing length. When the turbulence occurs in a thick layer we also have a sharp fluctuation of the mixing length. Increasing the mixing length increases the exchange coefficient dramatically. This feature has a close connection with the numerical stability of the scheme and with properties of positive inverse matrix operators.
This is the main problem for all methods. From results which are not completely presented, one can say that the best and most robust method is direct use of the nonlinear equation. This method still has limitations, but it can be improved.
The negative fluctuations effect of mixing length is increased by shallow convection. For this reason there is no possibility to solve the prognostic TKE equation with these methods when using the current shallow-convection scheme. One way to implement air saturation caused by turbulence is by adding the liquid water in the dependency in potential temperature. Without the shallow-convection scheme one can compute the tendency of TKE. This method does not use the complete TKE equation as a prognostic equation. This method is comparable to the current scheme for computation of fluxes, since it does not contain direct knowledge of "history" from the previous time-step. The difference is in the nonlinear approach of the computation.
For these properties we can solve at this moment TKE equation in "prognostic-diagnostic" regime with a very short integration step (around 10 s). To keep stability in appropriate range in all time-steps we initialize the TKE profile with the stationary TKE and then we compute the tendency of TKE with the described method. Some results are presented in the figures hereafter.
This method is good for scientific investigation, but not good for operational weather forecast. For longer time-step in case of rapid changes of the mean variables in the vertical direction, which can be caused by the discrete mesh (the density of mesh is insufficient), we obtain negative values for TKE below or above this region : which is not allowed. Without initialization with stationary TKE the model blows up after a few minutes.
To compare these methods for solving the nonlinear equation we did some experiments. These negative features can be suppressed by term-by-term splitting of the equation, by finding a better approximation for the derivatives or by changing the stability of the integration scheme. Other possibilities are to increase the number of vertical levels by additional interpolation. Interpolation can be done only for solving TKE. We have plan to investigate these ideas in the near future.
From analyses is apparent that the method for computation of fluxes depends on the distance to the surface. Kinematic turbulent fluxes are expressed close to the surface by TKE dependence. The computation assumes small change (10%) of fluxes in the ground layer. Therefore we apply constant statistical moments for this case. This approach is applied on the surface layer, which coincides with the last full level. In the free atmosphere we use a modified expression.
One can remark now, that spectral characteristics are used for the computation of some coefficients (cm, c e, ...), although the vertical grid is irregular. From linear spectral analysis, it is apparent that the integral value of total TKE is unaltered by a change of density of the grid. This method serves an estimation of the coefficients' magnitude. Their value, as was shown, must be tuned.
Every derivation is derived from the fact that TKE is computed on "half" levels. For the Lagrangian computation of horizontal advection it is more convenient to have TKE on "full" levels. For this reason we plan to compute the difference of TKE between two "half" levels.
For initialization the stationary (balanced) TKE is used.
Dry static energy for Richardson number expression is used in the ALADIN model. The moisture is included in the gas constant. The virtual potential temperature (there is a connection with dry static energy) is used for computation of the heat flux. Correction for shallow convection is not applicable at this moment.
For initialization of TKE we use a derived relation.
After taking into account our expressions and computing tendencies for TKE the exchange coefficients Kh , Km can be expressed. The computation of tendency for TKE is done in a diagnostic-prognostic regime at this moment. For the computation of the exchange coefficients we have to specify the parameterization coefficients c, c0 and c e. We suppose that cm is set from the spectral approach. Limitation of the Richardson number is garbled from current ALADIN scheme. The new form of the exchange coefficients, calculated at the time "+", repress the need for an anti-fibrillation scheme for the computation of mean variable tendencies.
During my 19 months stay in Bruxelles I improved my knowledge about Planetary-Boundary physics parameterization. I have a better overview of model ALADIN structure. I took part in the seminar and workshops. The complete scientific documentation of my 19 months stay in Brussels was created and distributed.