Solid Earth Open Access Open Access Open Access

Abstract. Meteorological numerical weather prediction (NWP) models solve a system of partial differential equations in time and space. Semi-lagrangian advection schemes allow for long time steps. These longer time steps can result in instabilities occurring in the model physics. A system of differential equations in which some solution components decay more rapidly than others is stiff. In this case it is stability rather than accuracy that restricts the time step. The vertical diffusion parametrization can cause fast non-meteorological oscillations around the slowly evolving true solution (fibrillations). These are treated with an anti-fibrillation scheme, but small oscillations remain in operational weather forecasts using ARPEGE and ALADIN models. In this paper, a simple test is designed to reveal if the formulation of particular a physical parametrization is a stiff problem or potentially numerically unstable in combination with any other part of the model. When the test is applied to a stable scheme, the solution remains stable. However, applying the test to a potentially unstable scheme yields a solution with fibrillations of substantial amplitude. The parametrizations of the NWP model ARPEGE were tested one by one to see which one may be the source of unstable model behaviour. The test identified the set of equations in the stratiform precipitation scheme (a diagnostic Kessler-type scheme) as a stiff problem, particularly the combination of terms arising due to the evaporation of snow.


Introduction
Meteorological numerical weather prediction (NWP) models solve a system of partial differential equations in time and space.Part of the atmospheric processes is resolved by the model dynamics, while the sub-grid, precipitation and radiation processes are said to be parametrized in the model physics.A number of years ago, eulerian advection schemes were replaced by the semi-lagrangian advection scheme in the model dynamics of the ARP ÉGE (Action de Recherche Petite Echelle Grande Echelle, which means research project on small-and large-scales) global model and ALADIN (Aire Limitée Adaptation Dynamique Développement International) limited area model.This has allowed the use of longer time steps and shortened the computer time needed to complete the operational weather forecast.Consequently these long time steps have been applied to the physics parametrizations, resulting in occasional instabilities in operational model forecasts using ARPEGE and ALADIN models.The physics tendencies are often larger than the dynamical ones, leading to strong transport in the vertical.
Previous studies (Kalnay and Kanamitsu, 1988;Girard and Delage, 1990;Bénard et al., 2000) have shown that the vertical diffusion parametrization can cause fast nonmeteorological oscillations, usually of 2 t period (where t is the model time step) around the slowly evolving true solution.These oscillations are referred to as fibrillations.The fibrillations are unwanted phenomena in the model since they might cause other model parts to produce false weather phenomena (clouds, rain, etc.) and lead to the wrong forecast or model blow-up in the operational suite.
The physical parametrization scheme in an operational NWP model often gradually evolves from one model version to the next one.Additional forcing terms are included in the particular parametrization, or the existing ones are made more complex to ensure better and more detailed representation of the real physical processes in the atmosphere.
Published by Copernicus Publications on behalf of the European Geosciences Union.
The numerical schemes used for parametrizations of physical processes did not achieve similar progress simultaneously (Beljaars, 1991).Numerically stable parametrization can destabilize when an additional forcing term is introduced (Beljaars, 2004).Consequently, more sophisticated parametrization can give worse results.
Instability can be the property of the equation system being solved, not the numerical scheme (Epperson, 2001).Such a system can be difficult to solve numerically even if it has bounded a continuous analytic solution for given initial and boundary conditions.Two strongly coupled processes of comparable timescales can result in solutions of different timescales and a stiff problem (Ropp et al., 2004).Atmospheric processes are strongly coupled so it is of little surprise that NWP models exhibit stiffness.
Stiff differential equations include several terms which vary at very different rates and can generate rapidly varying solutions.When explicit numerical methods for solving differential equations are applied to a stiff problem, an extremely small time step is required to maintain stability.Mathematica (Wolfram, 1996) uses an adaptive procedure to solve stiff differential equations, reducing the integration increment (time step) until a stable solution can be tracked.On the other hand, meteorological models sometime use a number of sub-steps in time to compute contributions of particular parametrization schemes (e.g.ECMWF, Beljaars, 2004).Stiffness is a transient phenomenon, creating difficulties for a mathematical criterion (e.g.Higham and Trefthen, 1993).
Stiffness is a property of the continuous system of equations, not of the numerical procedure used to solve it, but the instability that appears in its solution is referred to as numerical instability.These instabilities sometimes become evident as model "blow-up".Otherwise, the instability acts in a way that the model fields are noisy and oscillate in space and time.Although such oscillations are bounded and do not diverge too far from the true solution, they can be responsible for the false cloud or precipitation in the forecast (Teixeira, 2000).
One known stiff operator in meteorology is vertical diffusion.The non-linear diffusion equation can induce nonlinear instability through interaction of components of different wavelengths.The error is difficult to detect since the unstable solution only oscillates around the true one (Brown and Pandolfo, 1982).These oscillations are often called fibrillations and have a period of 2 t where t is the model time step.
The coupled system for vertical diffusion of momentum and heat was found unconditionally stable to perturbations for uniformly stratified and uniformly sheared flow (Davies, 1983).The source of instability in the heat diffusion equation is related to the formulation of the diffusion coefficient as a function of stability.
Linear stability analysis of the implicit scheme for vertical diffusion equation is unconditionally stable (Stull, 1988).However, the turbulent diffusion equations are nonlinear when the diffusion coefficient depends on the model variables.The practice of computing the exchange coefficient explicitly -from the current values of the model variables, and using this value for implicit computation of the future values of the model variables -leads to solutions that rapidly oscillate around a slow solution (Kalnay and Kanamitsu, 1988).The stability of the scheme decreases with increasing non-linearity of the diffusion coefficient.The solution can be stabilized using a filter that is applied only when the scheme is unstable.This filter is an over-implicit scheme that stabilizes the solution but reduces the accuracy.
The vertical diffusion schemes with non-linear diffusion coefficients are only conditionally stable.A solution has been proposed (Girard and Delage, 1990) to vary the implicitness of the scheme locally (in space and time).The scheme is over-implicit only where, when and as much as needed.Since the method used to stabilize the solution reduces the accuracy, it is applied only when the instability is detected through some criteria.The expressions for diffusion coefficients of heat and momentum can be different, yielding a more complex stabilization scheme.Bénard et al. (2000) implemented a stabilization scheme in the parametrization of turbulent diffusion in the ARP ÉGE and ALADIN models.Their stabilization scheme will be referred to as the antifibrillation scheme.The shallow convection scheme was also identified as a source of fibrillations.This was expected since it is based on the same principle as the turbulence scheme (Geleyn, 1987) and an associated stabilization scheme has been developed (Telišman et al., 1998).
An alternative solution to the problem of fibrillations is a construction of a numerical scheme designed for a particular problem.The 2 t oscillations were removed from the Met (Royal Meteorological Society) Office's UM (Unified Model) when a newly developed monotonically-damping second-order-accurate unconditionally-stable scheme was applied to the parametrization of vertical diffusion (Wood et al., 2007).The scheme also requires an estimate of the nonlinearity of the diffusion coefficient.Schlegel et al. (2012) implement a multi-rate time integration method to treat the stiff processes in the air pollution model.Teixeira (2000) studied the fibrillations and proposed another method for solving the diffusion equation.His method achieves stability by interpolating the data in the vertical and solving the equations on a new vertical grid that is determined by a prescribed stability criterion.There are other examples for numerical treatment of stiff equations from atmospheric chemistry (e.g.Jacobson, 2000).
Peer-reviewed literature seldom covers subjects such as failures of an operational forecast run or fibrillations.As will be shown here, there are 2 t oscillations of model fields that remain in the ARP ÉGE forecast although the anti-fibrillation schemes are activated in turbulence and shallow convection.Because these problems exist, they have motivated a search for a simple method that would detect a potentially unstable parametrization.The method is intended to test the model robustness and show if there is any parametrized process or a combination of processes that can produce numerically unstable results in any weather situation.The test should reveal a possibility of model failure or more moderate signs of numerical instability that could happen during the operational forecast.
In the absence of straightforward mathematical criteria for stiffness, the idea is to find a mechanism that would activate or amplify the numerically unstable behaviour if the parametrization being tested is prone to stiffness.The stiffness test of a potentially unstable scheme would result in model blow-up or significant amplification of fibrillations.Test of a stable scheme should yield a stable, non-fibrillating result.
The tests presented here were performed on the Météo-France ARP ÉGE model version that was operational until 2006.It is a global model that covers a wide variety of weather phenomena, and can use longer time steps.Both ARP ÉGE and ALADIN can use the same parametrization schemes.
The second section demonstrates how the idea of the stiffness test works on a simple reaction-diffusion equation.The third section describes its application in the ARP ÉGE model.The fourth section reveals which parametrized processes interact to produce amplified fibrillations.The last section presents a summary and conclusions.

Testing a simple non-linear diffusion equation 2.1 A simple non-linear diffusion equation
The stiffness test will first be applied to a simple non-linear damping equation as used by Kalnay and Kanamitsu (1988): where Kϕ p represents the non-linear diffusion coefficient with constant K = 10.The forcing term is D(t) = 1 − sin(2π t/24 h).The exponent p can be used to increase or decrease the non-linearity of the problem.Equation ( 1) is solved numerically using the method (g) from Kalnay and Kanamitsu (1988): (2) Method g from Kalnay and Kanamitsu (1988) is chosen for its similarity for the scheme used in ARP ÉGE and AL-ADIN.The diffusion coefficient Kϕ p is always computed explicitly.Parameter β is used to control the level of implicitness: for β = 0 the scheme is explicit, β = 0.5 yields a semiimplicit scheme, and for β = 1 the scheme is implicit.The IFS (Integrated Forecast System) model run operationally by ECMWF (European Centre for Medium-Range Weather   Forecasts) uses the over implicit scheme with β = 1.5.In ARP ÉGE, the anti-fibrillation scheme computes β for each grid point (Bénard et al., 2000).
Figure 1 shows solutions to Eq. ( 1) for p = 2 and t = 0.5 obtained with different levels of implicitness.The 2 t fibrillations are clearly visible in the solutions using β = 1 and β = 0.5, while the solution using β = 0 is unstable for this value of the time step.The amplitude of the fibrillations is an order of magnitude larger than the slowly varying referent solution.The fibrillations reduce in amplitude when the forcing term is small.Figure 2 shows solutions to the linear form of Eq. (1) (p = 0) and the non-linear problem (p = 2) with a shorter time step.For a linear problem (p = 0), implicit and semi-implicit schemes are stable and produce the correct solution (Fig. 2, left).When the time step used to solve the nonlinear form of Eq. (1) for p = 2 is reduced to t = 0.25, the implicit solution becomes free of fibrillations (Fig. 2 right).The explicit solution is always unstable, and the line escapes the area of the graph after a number of time steps.

The test results for the simple problem
In the stiffness test, the parametrization scheme that is subjected to the test is run with a different ( t/2) time step in the diffusion term (the first term on the right side of Eq. 1) than the one used in the "rest of the model" (other terms of Eq. 1).This way, the scheme produces a tendency of the model variable that will lead to inaccurate value of ϕ t+ t .That value will be used in the subsequent time step as input and produce a new tendency of the model variable, again using a modified time step in the diffusion term.The final solution is different from the reference one, and we are interested in its nature.
Figure 2 shows solutions to a linear form of Eq. ( 1   non-linear form with p = 2 on the right side.The results of the proposed test applied to these two problems are shown in Fig. 3.When the test is applied to a linear problem, the result remains non-oscillatory with a shift in amplitude of the slowly varying solution for the implicit and semiimplicit schemes, while the explicit scheme is always unstable (Fig. 3, left).Applying the test to a non-linear problem yields a solution where fibrillations dominate over the slowly evolving component (Fig. 3, right).The solution that behaved properly in the reference run (the implicit one) exhibits strong fibrillations with amplitude an order of magnitude larger than the true solution (Fig. 3).
It is shown here that when the test is applied to a stable scheme, the solution remains stable (although offset from the true solution), but applying the test to a potentially unstable scheme yields a solution with fibrillations of substantial amplitude.

Model description
ARP ÉGE is a global spectral model (Courtier and Geleyn, 1988) used for operational forecast by Météo-France (Courtier et al., 1991).Its limited area version, called AL-ADIN, is run operationally by many national weather services, including the Croatian one (Ivatek-Šahdan and Tudor, 2004).The model version used here is a semi-implicit semi-lagrangian two-time-level scheme (Robert, 1982) for the hydrostatic shallow atmosphere.The non-linear residual is computed using the stable extrapolation for the semilagrangian scheme (Hortal, 2002).The tendencies of the model variables due to parametrized processes are computed before the dynamical ones and interpolated to the origin point of the semi-lagrangian trajectory.Both models use hybrid η coordinate (Simmons and Burridge, 1981) and finite differences in the vertical.ALADIN uses double Fourier representation of the fields (Machenhauer and Haugen, 1987).
The computations of the parametrized processes are done for a single vertical column and compute the vertical fluxes of model variables.The vertical flux of variable ψ due to the parametrized process is F proc ψ .It is computed as a function of model variables, their vertical gradients, other parameters and the model time step: (3) The fluxes and the model forecast variables are computed on a grid that is staggered in the vertical.The model variables are computed on model levels, while the fluxes are computed on intermediate levels between the model levels.The model grid is not staggered in the horizontal.All the parametrization schemes use the same input values of the model variables, the computations are done in parallel and the physics tendency is obtained from the fluxes due to different processes at the end of the physics computations.The local physical tendencies are computed as vertical divergence of the fluxes.The partial tendency of variable ψ due to effects of parametrized physical processes is where the square brackets contain all parametrized processes relevant for the evolution of the model variable ψ (e.g.turb   for vertical turbulent diffusion, prec for resolved precipitation, conv for deep convection, etc.).The horizontal effects are considered only through moisture convergence used in the deep convection scheme.Other horizontal effects -the exchange between columns -is treated by the dynamical part of the model through advection and horizontal diffusion.
The model version used in this study is based on the following hypotheses of atmospheric energy and conservation aspects: the atmosphere is composed of a mixture of two perfect gasses (dry air and water vapour), and the system is in thermodynamic equilibrium (Courtier and Geleyn, 1988).In the tested version of the model, the condensed phases that appear in the atmosphere have zero volume, fall and exchange temperature and momentum with the layers they cross, and leave the system in one time step, removing mass, momentum and energy.
Temperature and specific moisture are prognostic variables in the model, but parametrizations compute fluxes of heat.The temperature tendency is computed from heat fluxes using specific heat at constant pressure, gas constant and latent heats of evaporation and sublimation.Specific heat at constant pressure c p varies with specific moisture q v , but remains independent of temperature: where c pa is the specific heat for dry air and c pv is the specific heat for water vapour.The gas constant R varies with specific moisture q v as where R d is the gas constant for the dry air and R v is the gas constant for water vapour.Evaporation and sublimation latent heats vary with temperature: where L 0 is the latent heat at the triple point temperature T t , c l is the specific heat of liquid water and c i is the specific heat of ice.

Stiffness tests of the parametrization schemes
The model parametrizations are designed to restore the thermodynamic equilibrium during one time step using the parametrized processes.In the test, the flux is computed using the modified time step, e.g.
The original model time step is used in subsequent computation of the model tendency due to parametrized processes.The computed flux is modified and the model does not reach the desired balanced state exactly.
As a consequence, the model might drift from the true solution, and the forecast will be incorrect, producing stronger or weaker temperature diurnal cycle, more or less cloudiness and rain and more or less momentum transport in the vertical.Any of these effects are an expected consequence of the imbalance imposed by the modified time step used in the particular flux computation.We are interested to see if the model variables exhibit unstable behaviour during the forecast run, such as fibrillations.
The appearance of fibrillations in the model fields is diagnosed using the amplitude of the 2 t oscillations that is computed on the model levels using values of the model variables from three subsequent time steps.The magnitude of fibrillations in the temperature field is computed using the formula For simplicity, only the results for temperature are presented here.The ARP ÉGE model is run with a time step of 830.77 s, equivalent to 13 time steps per 3 h.The model used a stretched grid of variable horizontal resolution.The resolution is highest above France and lowest on the opposite side of the globe.In the vertical the model employs 41 terrains following η levels; the lowest level lies approximately 17 m above the surface.The lowest level is numbered 41 and the level numbers decrease with height.
The model is run for a 96 h period, with the magnitude of fibrillations computed according to Eq. ( 9).A map of the amplitude of fibrillations of temperature from the reference model run (without the test) after 96 h forecast on the lowest model level is shown in Fig. 4. The areas with amplitude larger than 0.5 Kelvin are shaded.Fibrillations exceed the 0.5 K threshold in a number of areas.Figure 5 shows the evolution of temperature on the 12 lowest model levels for each time step in one model grid point (lon = 7.22 • E, lat = 60.02• N) during the 96 h forecast of the reference model run.The fibrillations are transient in nature and active only in certain weather conditions.In this model grid point, the fibrillations occur in stable conditions with a temperature inversion.
The test was performed for each parametrization scheme, one by one.In each test, the model was run with the operational time step, and the time step used in the tested parametrization was modified.This way, the parametrization computed the vertical flux of the model variable that is needed to "balance" the model state in half the time step.This flux is combined with fluxes from other parametrizations (computed with the operational time step) and the resulting tendency of the model variable and the future state of the model at t + t are all computed using the operational time step.As noticed by one of the reviewers, any perturbation is quickly damped for an unconditionally stable numerical scheme, with strong damping implying that if a good  The amplitude of the oscillations of the model variable from the reference model run (Figs. 4 and 5) is compared to the amplitudes from the run with the test (for example Figs. 6 and 7) to compare the noise when the error is introduced into the parametrization being tested.Figure 6 shows the amplitude of 2 t oscillations of temperature computed according to Eq. ( 9) on the lowest model level after 96 h period for the model run, where the test with the amplification of fibrillations was activated in the large-scale precipitation scheme.The amplitude exceeds 16 K in several areas.One should keep in mind that ARP ÉGE uses stretched grid with highest resolution above France and an order of magnitude lower resolution on the farther side of the globe.Figure 4 shows the amplitude plotted as a shaded point for each model grid point.Consequently, the area of considerable amplification in the amplitude close to North Pole, north of the Bering Sea, appears as stripes of grid points, although this is a continuous area in the model.Figure 7 shows the evolution of temperature on the lowest model level for the same test in the large-scale precipitation scheme.The amplitude of fibrillations varies in time for several orders of magnitude in a particular model grid point.The fibrillations appear, amplify and diminish during the model forecast, and the model continues to run even though the variations in the temperature field are far from possible to be achieved in nature.
The test results for each parametrization scheme can be summarized as follows: -Vertical turbulent diffusion and shallow convection are already known to be stiff.These parametrizations already have the anti-fibrillation scheme implemented.19 Fig. 6.The amplitude of t temperature oscillations on the lowest model level (in Kelvin) after 96 h forecast, computed according to Eq. ( 9) for the model run with the fibrillation amplification test of the large-scale precipitation scheme.
-Tests of the deep convection and gravity wave drag schemes did not produce significant changes in the magnitude of fibrillations.
-The test of the large-scale precipitation scheme produced high amplitude fibrillations in some parts of the model domain (Fig. 6).The amplitude of fibrillations (Fig. 7) increased several orders of magnitude.
The test results revealed another stiff process in the model -the large-scale precipitation scheme.The scheme describes processes of condensation, evaporation, freezing and melting of precipitation.It also interacts with other processes described in the model.It was important to find which parametrized process, or a combination of processes, was responsible for the amplification of fibrillations.
4 The source of instability in the stratiform precipitation scheme

The stratiform precipitation scheme
The stratiform precipitation scheme computes the precipitation fluxes due to resolved condensation.The simplified Kessler-type scheme (Kessler, 1969)   of a grid cell.The specific water vapour of the saturated air at the wet bulb temperature q w has a symmetric role, being the lower limit for the specific humidity q v in the case of condensation and the upper limit in the case of evaporation.Saturation is treated using the difference between the wet bulb and local grid-average moisture (q w − q v ): -if q w −q v < 0 there is condensation in the layer and only the amount q = q v − q w is condensed, the condensate leaves the layer as precipitation flux and the layer becomes saturated, -if q w − q v > 0 and the precipitation flux on the top of the layer is P t > 0, the precipitation evaporates, but not more than is needed to saturate the layer.
The precipitation flux P changes in the vertical are computed using the pressure coordinate p: where the evaporation of precipitation is computed as C evap = 4.8 × 10 6 is the value of the tunable parameter, r m is the snow to water ratio, and R evgsl = 80 is the ratio of evaporation speeds for snow and rain.The assumptions used in the large-scale precipitation scheme are the following: -the precipitation is at the same temperature as the surrounding air, -it is all liquid if the air temperature is above the triple point temperature T t and all frozen if it is below T t .In the case of condensation, we must consider if the temperature is warmer or colder than the triple point: -if T > T t , the excess water vapour condenses as liquid and the snow to water ratio of the precipitation flux decreases by r m = r m P t P b where P t and P b are the precipitation fluxes at the top and at the bottom of the layer.
-if T < T t , the excess water vapour condenses as solid and the snow to water ratio of the precipitation flux increases r m = 1 − (1 − r m ) P t P b .Melting (T > T t ) and freezing (T < T t ) of the precipitations are computed as the change of the snow to water ratio in the vertical: where the melting and freezing of precipitations is computed as and C melt = 2.4 × 10 4 is the value of the tunable parameter.
The final precipitation fluxes of liquid and solid condensates are computed according to where P L is the precipitation flux for the liquid and P S for the solid condensates.

The search for the source of the instability
The next issue is to find which process in the model, or a combination of processes, is responsible for amplification of fibrillations in the stiffness test of the large-scale precipitation scheme.Here it is assumed that fibrillations can be removed or reduced in intensity if one of the processes responsible for amplification is switched off or modified in intensity.For this purpose, a number of experiments are executed in which the model is run with the stiffness test in the largescale precipitation scheme, and one of the parametrization schemes is switched off (or reduced in intensity).The results of the experiments that test the interactions with other parametrization schemes can be summarized as follows: -the shallow convection scheme is switched off, but fibrillations remain, -the deep convection is switched off, and the fibrillations amplify, especially above sea surface on winter hemisphere, -switching off the vertical diffusion scheme leads to model blow-up, -the gravity wave drag is switched off, but the fibrillations remain, -the radiation scheme is switched off, but the fibrillations remain.
These tests show that the fibrillations are not the result of coupling a stratiform precipitation scheme to other parametrization schemes (even those that are known to be stiff).
The large-scale precipitation scheme includes several processes.These processes were modified by means of tuning parameters in separate experiments.In each experiment, one process was modified, and the other processes remained unchanged.Figures 8 to 13 show maps of amplitude of  9) for the model run with the fibrillation amplification test of the large-scale precipitation scheme but with reduced freezing of precipitation.The figure on the right shows temperature (in Kelvin) in one model point on the 12 lowest model levels for each time step during the 96 h forecast for the model run with the fibrillation amplification test of the large-scale precipitation scheme but with reduced freezing of precipitation.fibrillations of temperature on the lowest model level after 96 h integration, and the evolution of temperature on the 12 lowest model levels in one grid point for the experiments that examine the processes in the large-scale precipitation scheme.The results of these experiments are the following (related figures are listed in the brackets): -the coefficient for evaporation of precipitations (C evap ) is set to 0, and the fibrillations disappear (Fig. 8), -the coefficient for melting of precipitations (C melt ) is significantly reduced, and the fibrillations remain (Fig. 9), -the condensation is turned off through logical switch, and the fibrillations disappear (Fig. 10), -the cryoscopic cycle is turned off through logical switch, and the fibrillations disappear (Fig. 11), -the cryoscopic cycle is activated, but the coefficients for evaporation (C evap ) and melting (C melt ) of precipitations are set to 0, and the fibrillations disappear completely.
Fibrillations are linked to evaporation and the cryoscopic cycle.One parameter that affects both processes is the ratio of the speed of evaporation of precipitation between ice and water (R evgsl ) that is operationally set to 80.
Large values of the ratio of the speed of evaporation of ice and water are the primary cause of the fibrillation amplification acting only when the cryoscopic cycle is activated together with the evaporation of precipitation.9) for the model run with the fibrillation amplification test of the large-scale precipitation scheme but with evaporation ratio of ice and water reduced to 8. The figure on the right shows temperature (in Kelvin) in one model point on the 12 lowest model levels for each time step during the 96 h forecast for the model run with the fibrillation amplification test of the large-scale precipitation scheme but with evaporation ratio of ice and water reduced to 8.
The large value of the coefficient describing the ratio of the speed of evaporation of ice and water is a consequence of several assumptions: -The first assumption states that all the condensed water and ice falls to the ground during one model time step.
-The second assumption is that the condensates are of the same temperature as the surrounding air, and it is liquid if the temperature is above the triple point, and solid if it is below.
-Water drops fall about 4 times faster than ice crystals of the same weight (Pielke, 2002).Both fall through all the layers below the cloud during one time step.It is therefore assumed that the ice crystals spend 4 times more time in a certain layer.Slower fall speed is accounted for by adapting the parametrization so that the ice crystals evaporate 4 times faster than the water drops.
-Additionally, one should take into account that ice crystals have much larger surface than water droplets of the same mass, so they can evaporate much faster.This is accounted for by increasing the ratio of evaporation rates for ice and water even further.
When all factors are considered, the ratio of the speed of evaporation of ice crystals and water droplets is set to 80. -setting R evgsl = 1 means that snow evaporates as rapidly as rain -there are no fibrillations, -setting R evgsl = 20 leads to small reduction of fibrillation magnitude, -setting R evgsl = 8 removes fibrillations from certain areas (Fig. 12), -setting R evgsl = 4 removes fibrillations from most areas (Fig. 13).
When the ice crystal falls through the unsaturated warm layer, it evaporates and melts rather quickly by using the heat from the layer and reducing the temperature of the layer.In the next time step, the layer becomes colder and saturated so the ice crystals that fall through the layer grow, removing water vapour from the layer.Resulting condensation produces latent heating and increases its temperature.In the next time step, the layer is unsaturated and warm again.
The radiation scheme computes new fluxes every time step so it could interact with the instabilities caused by the largescale precipitation scheme if there are alternating time steps with and without a cloud.However, the radiation contribution to the evolution of temperature during one time step is small when compared to contributions of other parametrizations, such as large-scale precipitation or vertical diffusion.This is why it does not have a significant affect on the instability produced by the large-scale precipitation scheme.

Conclusions
Noisy solutions are a common feature of atmospheric models.Girard and Delage (1990) apply the stabilization scheme for vertical diffusion in the Canadian spectral weather forecast model, Bénard et al. (2000) in ARP ÉGE and ALADIN models, and Wood et al. (2007) develop a stable scheme for the Unified Model.Teixeira (2000) describes noisy solutions in the IFS model that are a consequence of the interaction between the radiation, vertical diffusion and cloud schemes.
In a meteorological model, the equations that parametrize physical processes are highly non-linear and depend on model prognostic variables and their derivatives.Several parametrizations describe processes that have timescales comparable or even less than the model time step.The problem of non-linearity is made worse in parametrizations that use different formulations or coefficients depending on the weather conditions, such as the static stability of the boundary layer or the condition if the air temperature is above or below the triple point.
Here it is shown how to test the parametrizations for stiffness as they are implemented in the model.The parametrizations are "pushed" to become unstable by the modified time step used only in the particular computations.The test activates and amplifies the fibrillations.When the test was applied to a non-stiff problem, no fibrillations appeared.But when the test was applied to a stiff system, the fibrillations were initiated and amplified by several orders of magnitude.In meteorological model, the fibrillations amplified only in certain areas of the globe.The fibrillations were also transient in time, confirming that the fibrillations are activated only when particular conditions are met.
The method presented here provides a technical test of the model robustness that checks for instability resulting from stiffness.Using the presented test revealed another stiff equation set in the model, the one describing resolved precipitation processes.The findings of this paper support the previous results of Brown and Árnason (1973) and Brown (1977) who showed the equations of droplet growth to be stiff.Similarly, Sednev and Menon (2012) analysed source codes of bulk cloud microphysics parametrizations for warm rain processes in various community models and found deficiencies in numerics that restrict the length of the time step.These equation sets describe growth of droplets of different sizes, and are more complex than the simple parametrization scheme used here.
Explicit numerical methods are not useful when solving stiff sets of equations.The equations for droplet growth can be solved using an iterative (Brown and Árnason, 1973) or a non-iterative partially implicit solution (Brown, 1977).Other examples of the numerical methods for treatment of stiff equations are those used for the turbulent diffusion and to describe chemical properties (Jacobson, 2000).

Fig. 1 .
Fig. 1.Numerical solutions of Eq.(1 with implicit (dotted), trapezoid (short dash line) and explicit (long dash line) schemes for exponent p = 2 and time step t = 0.5 h.Function D(t) is also shown (solid line) for reference.

Fig. 2 .
Fig. 2. The left figure shows the referent numerical solution of the equation 1 (using scheme described by equation 2) for the linear problem with p = 0.The right figure shows the solution to the non-linear problem with p = 2 and time-step ∆t = 1/4 hours.The results are shown for implicit (dotted), trapezoid (short dash) and explicit (long dash) schemes.Function D(t) is also shown (full line) for reference. 17

Fig. 2 .
Fig. 2. The left figure shows the referent numerical solution of Eq. (1) (using scheme described by Eq. 2) for the linear problem with p = 0.The right figure shows the solution to the non-linear problem with p = 2 and time step t = 1/4 h.The results are shown for implicit (dotted line), trapezoid (short dash line) and explicit (long dash line) schemes.Function D(t) is also shown (solid line) for reference.

Fig. 3 .
Fig. 3.The left figure shows the result of solving the equation 1 (using scheme described by equation 2) with the fibrillation amplification test for the linear problem with p = 0 and the right figure shows the solution to the non-linear problem with p = 2 and time-step ∆t = 1/4 hours.The results are shown for implicit (dotted), trapezoid (short dash) and explicit (long dash) schemes.Function D(t) is also shown (full line) for reference.

Fig. 4 .
Fig. 4. The amplitude of 2∆t temperature oscillations on the lowest model level (in Kelvin) after 96 hour forecast for the referent model run computed according to Equation 9.

Fig. 3 .
Fig. 3.The left figure shows the result of solving Eq. (1) (using scheme described by Eq. 2) with the fibrillation amplification test for the linear problem with p = 0, and the right figure shows the solution to the non-linear problem with p = 2 and time step t = 1/4 h.The results are shown for implicit (dotted line), trapezoid (short dash line) and explicit (long dash line) schemes.Function D(t) is also shown (solid line) for reference.

Fig. 4 .Fig. 4 .
Fig. 4. The amplitude of 2∆t temperature oscillations on the lowest model level (in Kelvin) after 96 hour forecast for the referent model run computed according to Equation 9.

Fig. 5 .Fig. 5 .
Fig. 5. Temperature (in Kelvin) in one model point on the 12 lowest model le 96 hour forecast for the referent model run.Level 41 is the lowest model leve

Fig. 6 .
Fig. 6.The amplitude of 2∆t temperature oscillations on the lowest model level (in Kelvin) after 96 hour forecast computed according to Equation 9 for the model run with the fibrillation amplification test of the large scale precipitation scheme.

Fig. 7 .Fig. 8 .Fig. 7 .
Fig. 7. Temperature (in Kelvin) in one model point on the 12 lowest model lev 96 hour forecast for the model run with the fibrillation amplification test of the

Fig. 8 .Fig. 8 .
Fig. 8.The left figure shows the amplitude of 2∆t temperature oscillations on the lowest model level (in Kelvin) after 96 hour forecast computed according to Equation 9 for the model run with the fibrillation amplification test of the large scale precipitation scheme but when the evaporation of precipitation is switched off.The figure on the right shows temperature (in Kelvin) in one model point on the 12 lowest model levels for each time-step during the 96 hour forecast for the model run with the fibrillation amplification test of the large scale precipitation scheme but when the evaporation of precipitation is switched off.

Fig. 9 .Fig. 10 .Fig. 9 .
Fig. 9.The left figure shows the amplitude of 2∆t temperature oscillations on the lowest model level (in Kelvin) after 96 hour forecast computed according to Equation 9 for the model run with the fibrillation amplification test of the large scale precipitation scheme but with reduced freezing of precipitation.The figure on the right shows temperature (in Kelvin) in one model point on the 12 lowest model levels for each time-step during the 96 hour forecast for the model run with the fibrillation amplification test of the large scale precipitation scheme but with reduced freezing of precipitation.

Fig. 9 .Fig. 10 .Fig. 10 .
Fig. 9.The left figure shows the amplitude of 2∆t temperature oscillations on the lowest model level (in Kelvin) after 96 hour forecast computed according to Equation 9 for the model run with the fibrillation amplification test of the large scale precipitation scheme but with reduced freezing of precipitation.The figure on the right shows temperature (in Kelvin) in one model point on the 12 lowest model levels for each time-step during the 96 hour forecast for the model run with the fibrillation amplification test of the large scale precipitation scheme but with reduced freezing of precipitation.

Fig. 11 .Fig. 12 .Fig. 11 .
Fig. 11.The left figure shows the amplitude of 2∆t temperature oscillations on the lowest model level (in Kelvin) after 96 hour forecast computed according to Equation 9 for the model run with the fibrillation amplification test of the large scale precipitation scheme but with cryoscopic cycle switched off.The figure on the right shows temperature (in Kelvin) in one model point on the 12 lowest model levels for each time-step during the 96 hour forecast for the model run with the fibrillation amplification test of the large scale precipitation scheme but with cryoscopic cycle switched off.

Fig. 11 .Fig. 12 .Fig. 12 .
Fig. 11.The left figure shows the amplitude of 2∆t temperature oscillations on the lowest model level (in Kelvin) after 96 hour forecast computed according to Equation 9 for the model run with the fibrillation amplification test of the large scale precipitation scheme but with cryoscopic cycle switched off.The figure on the right shows temperature (in Kelvin) in one model point on the 12 lowest model levels for each time-step during the 96 hour forecast for the model run with the fibrillation amplification test of the large scale precipitation scheme but with cryoscopic cycle switched off.

913, 2013 908 M. Tudor: Numerical instability and stiffness of the ARP ÉGE and ALADIN models
www.geosci-model-dev.net/6/901/2013/Geosci.Model Dev., 6, 901- In the last set of experiments, the original value of R evgsl is changed from 80 to other values and the results are as follows: The left figure shows the amplitude of 2∆t temperature oscillations on the lowest model level (in Kelvin) after 96 hour forecast computed according to Equation 9 for the model run with the fibrillation amplification test of the large scale precipitation scheme but with evaporation ratio of ice and water reduced to 4. The figure on the right shows temperature (in Kelvin) in one model point on the 12 lowest model levels for each time-step during the 96 hour forecast for the model run with the fibrillation amplification test of the large scale precipitation scheme but with evaporation ratio of ice and water reduced to 4. The left figure shows the amplitude of 2 t temperature oscillations on the lowest model level (in Kelvin) after 96 h forecast computed according to Eq. (9) for the model run with the fibrillation amplification test of the large-scale precipitation scheme but with evaporation ratio of ice and water reduced to 4. The figure on the right shows temperature (in Kelvin) in one model point on the 12 lowest model levels for each time step during the 96 h forecast for the model run with the fibrillation amplification test of the large-scale precipitation scheme but with evaporation ratio of ice and water reduced to 4.