Generalization and application of the flux-conservative thermodynamic equations in the AROME model of the ALADIN system

General yet compact equations are presented to express the thermodynamic impact of physical parameterizations in a NWP or climate model. By expressing the equations in a flux-conservative formulation, the conservation of mass and energy by the physics parameterizations is a built-in feature of the system. Moreover, the centralization of all thermodynamic calculations guarantees a consistent thermodynamical treatment of the different processes. The generality of this physicsdynamics interface is illustrated by applying it in the AROME NWP model. The physics-dynamics interface of this model 5 currently makes some approximations, which typically consist of neglecting some terms in the total energy budget, such as the transport of heat by falling precipitation, or the effect of diffusive moisture transport. Although these terms are usually quite small, omitting them from the energy budget breaks the constraint of energy conservation. The presented set of equations allows to get rid of these approximations, in order to arrive at a consistent and energy-conservative model. A verification in an operational setting shows that the impact on monthly-averaged, domain-wide meteorological scores is quite neutral. However, 10 under specific circumstances, the supposedly small terms may turn out not to be entirely negligible. A detailed study of a case with heavy precipitation shows that the heat transport by precipitation contributes to the formation of a region of relatively cold air near the surface, the so-called cold pool. Given the importance of this cold pool mechanism in the life-cycle of convective events, it is advisable not to neglect phenomena that may enhance it.


Introduction
The conservation of mass and energy are important characteristics of a numerical atmospheric model.Especially in view of the application in climate studies, even small violations of the conservation laws can accumulate over a long integration time, and lead to faulty results (Staniforth and Wood, 2008;Lucarini and Ragone, 2011).Atmospheric forecast models are usually constructed by combining a dynamical core with physical parameterizations.In general, the dynamical core describes the atmospheric behaviour up until the resolved scales, while the physical parameterizations estimate the effect of subgrid processes (Gassmann, 2013).
A lot of research has been spent in designing dynamical cores that conserve mass and energy (Thuburn, 2008).Common strategies include a careful selection of the prognostic variables (Ooyama, 1990(Ooyama, , 2001;;Klemp et al., 2007), the formulation of the equations in flux form (Satoh, 2003), or taking advantage of properties of the Hamiltonian character of the atmospheric equations (Salmon, 2004;Gassmann and Herzog, 2008;Zängl et al., 2015).In contrast with these efforts on the dynamical core, the energy conservation and consistent thermodynamics seem to be less of a priority in the development of the physical parameterizations, or in the way they are coupled to the dynamical core.
A possible explanation is that the thermodynamics of the dynamical core are less complicated than those of the physical parameterizations.More specifically, the dynamics are usually considered adiabatic and reversible (except for numerical diffusion) (Gassmann, 2013).The physics parame-D.Degrauwe et al.: Generalization and application of the flux-conservative thermodynamic equations terizations, on the other hand, include mass and energy exchange with the surface, as well as radiative fluxes at the top of the atmosphere.They constitute an open thermodynamic system, for which the conservation laws are more difficult to enforce.Moreover, it is tempting to consider physics parameterizations as plug-compatible, i.e. they are considered as a black box which, given an atmospheric state, returns an effect on the dynamical prognostic variables.Unfortunately, this plug compatibility seems to go at the expense of carefully investigating the thermodynamic consistency between the dynamical core and the physics parameterizations, and inserting a new parameterization in a model comes with implicit assumptions and ad-hoc approximations.
There is, however, an increased interest in different aspects of the coupling of physical parameterizations to the dynamical core.One of the issues is the organization of the time step.This problem has been studied with academic toymodels (see, e.g.Caya et al., 1998;Staniforth et al., 2002;Termonia and Hamdi, 2007), as well as in 3-D models (Hortal, 2002;Williamson, 2002).The thermodynamic aspects of the physics-dynamics coupling is another topic that deserves some attention.Although some attempts have been made to rigorously formulate the equations for a multicomponent atmosphere (Ooyama, 2001;Bannon, 2002), it remains a fact that many operational models make several adhoc approximations (Bryan and Fritsch, 2002).Catry et al. (2007), hereafter CGTBT07, presented a set of equations that expresses the effects of physics parameterizations in a fluxconservative formulation.The advantage of this approach is that this is an inherently mass-and energy-conservative system.
The current paper develops the proposal of CGTBT07 further by generalizing it for a system with an arbitrary number of hydrometeors with arbitrary interactions between them.It should be emphasized that the scope of this work is limited to the coupling of the atmospheric physics parameterizations to the dynamical core.For instance, when energy-conserving equations are presented, this property does not necessarily hold for the atmospheric model as a whole, but only regarding the influence of the physical parameterizations.Other aspects of the model, most notably its dynamical core, may not be energy conserving.Also the mutual interactions between different parameterizations are not considered in this paper, as they relate only indirectly to the time evolution of the prognostic atmospheric variables.The next section presents the equations of this generalized system.In Sect.3, this set of equations is applied in the AROME numerical weather prediction (NWP) model (Seity et al., 2011), thus allowing to get rid of some approximations that are currently made.Section 4 discusses the impact on the meteorological results, both by means of monthly scores and with an in-depth case study of a cold pool formation under heavy precipitation.Section 5 presents the conclusions.
2 Formulation of the generalized flux-conservative equations

Framework of hypotheses
Because the behaviour of the atmosphere is too complex to be described exactly, every numerical model needs to make simplifying hypotheses.This is no different for the work described in the current paper.It is not our aim to present a set of equations which is exact in the sense that it is free of approximations.However, a crucial aspect of the work presented in CGTBT07 is that the set of hypotheses that relate to the thermodynamics is defined from the very beginning.This is important for two reasons.First, it ensures that the simplifications act consistently throughout the model.Second, it allows to set some non-negotiable constraints.For instance, the conservation of energy must be satisfied, no matter what other simplifications are made.This approach of setting the simplifying hypotheses from the beginning contrasts with the conventional approach of ignoring supposedly small terms along the way.
The framework of hypotheses is the following: -A fully barycentric view of air parcels is adopted.This means that all hydrometeors (both suspended and precipitating) are considered as integral parts of the air, and contribute to the parcel's motion, density, and heat capacity.This barycentric view has been studied and motivated by many researchers (Wacker and Herbert, 2003;Bott, 2008;Gassmann and Herzog, 2008).
-Water condensates are assumed to have zero volume.This is a common approximation in atmospheric modelling.
-Temperature is homogeneous across all species, even falling hydrometeors.For small hydrometeors, this approximation is easily justified, given their short relaxation time (Bott, 2008).For larger hydrometeors, it is a cruder approximation, but it goes together with the barycentric view: since such hydrometeors are considered part of the parcel, they also take the parcel's temperature.
-The specific heat values of all species are constant with temperature.
-The latent heat values of sublimation and evaporation, L i and L l , respectively, vary linearly with temperature T : where T 0 = 0 K, c pv is the specific heat capacity at constant pressure of water vapour, and c i and c l are the specific heat capacity values of ice and liquid water, respectively.
It should be mentioned that this same framework of assumptions has been used by Marquet (2011), Marquet andGeleyn (2013), andMarquet (2015) to cleanly develop moist atmospheric thermodynamic quantities such as moist entropy, moist potential temperature, and moist Brunt-Väisälä frequency.

2.2
The flux-conservative equations for a system with five water species The system considered in CGTBT07 consists of dry air (specific mass fraction q d = ρ d /ρ tot ) plus five prognostic water species: vapour (specific mass fraction q v ), suspended liquid water droplets (q l ), suspended ice crystals (q i ), precipitating rain (q r ), and precipitating snow (q s ).For this system, the following equations are derived for the time evolution of the prognostic species due to physical parameterizations: In these equations, P k denotes precipitation fluxes and D k denotes diffusive fluxes.Note that it is necessary that k=d,v,i,l D k = 0 to ensure that all terms on the right-hand sides cancel out.The terms R k 1 ,k 2 denote pseudofluxes and represent mass transfer between two water species.The concept of pseudofluxes is essential to the presented system and deserves some more explanation.The common and more intuitive way to express a mass transfer between two species is through a time tendency.For instance, consider the microphysical process of condensation, which is a mass transfer from water vapour to liquid cloud water droplets.The effect of this process on the specific humidities could be expressed as The pseudoflux R v,l expresses exactly the same effect, only as a flux instead of as a tendency.This flux is determined by taking the vertical integral of the tendency: Although a pseudoflux is arguably more difficult to interpret than a tendency, writing conversions between species in terms of pseudofluxes offers the possibility to write the evolution equations in a flux-conservative form.The benefit of this is explained further.Also note that this does not mean that the internals of the physics parameterizations should be formulated in terms of pseudofluxes.Instead, it is only at the moment when the contributions of the physics parameterizations are added to the prognostic variables, that pseudofluxes are beneficial.They can be determined at that point from the more conventional tendencies using the expression above.
The thermodynamic equation for the system with 4 hydrometeors is as follows: where ĉ = c pd q d +c pv q v +c i q i +c l q l 1−q r −q s , and J s and J rad are the diffusive and radiative heat fluxes, respectively.c p is the total heat capacity of the parcel, given by c p = c pd q d + c pv q v + c l (q l + q r ) + c i (q i + q s ).
It should be noted that Eq. ( 8) expresses only the thermodynamic effect of the physical parameterizations.The complete thermodynamic equation of the atmospheric model would also include terms that are resolved by the dynamics of the model.
A full discussion of these equations is given in CGTBT07, but we would like to stress the following characteristics: -All equations are flux-conservative, i.e. every right-hand side is a divergence of a summation of fluxes.The importance of this property cannot be underestimated, because it means that this system intrinsically conserves mass and energy.Put somewhat simplistically, in a fluxconservative system, the only way energy or mass can leave one model layer, is by transporting it to an adjacent layer.Therefore, mass and energy are conserved by design of the system.
To derive these relations, one starts from the definition of a flux as a product of a density with a velocity.For instance, for rain, one writes The absolute velocity w * of the center of mass of the parcel is given by the weighted average of the velocities of the components.In a system where only rain and snow are precipitating, this means that The relative velocity of rain is then given by w r = w * r − w * , so the relative precipitation flux becomes -The latent heat values of sublimation and condensation L i and L l that appear on the right-hand side, are evaluated at T 0 = 0 K.This does not mean that the temperature dependency of these latent heat values is neglected.Instead, it is accounted for by considering the time derivative of the enthalpy c p T .Considering only the process of condensation, the traditional way to express its thermodynamic effect would be Using the before-set assumption that L l varies linearly with temperature, and the fact that, still only considering condensation, ∂q v /∂t = −∂q l /∂t, so ∂c p /∂t = (c l − c pv )∂q l /∂t, this expression becomes which can be rewritten as This shows how the temperature dependence of the latent heat values can be accounted for by considering the tendency of enthalpy.
-Although the equations only describe the evolution of water species, similar flux-conservative equations could be formulated for other atmospheric variables like momentum, turbulent kinetic energy, etc.In this paper, only water species and their effect on the thermodynamic equation are studied.

The generalized flux-conservative equations
Despite the clear strength of the equations proposed by CGTBT07, their application is not straightforward because of the fixed number of water species, and because of the fixed set of interactions between them (six pseudofluxes).More advanced microphysics schemes often consider more water species, for instance by including graupel and/or hail (Lascaux et al., 2006), or by separating convective and nonconvective fractions of hydrometeors (Piriou et al., 2007).Also the fact that only six transfer mechanisms between the water species are possible is limiting.For instance, snow melting cannot be represented directly, but it should be written as a combination of snow sublimation (R s,v ) and rain evaporation (R r,v ).Although thermodynamically fully correct, it would be better to have a system that digests all kinds of transfers between water species.It is, however, possible to generalize the equations from CGTBT07, without touching the important characteristics.We introduce the following notation: n is the number of water species, the index k = 1, . .., n denotes a single water species, and by convention, k = 0 denotes the dry air component.The specific heat capacity values at constant pressure of the different species are written generically as c k , the latent heat of evaporation or sublimation at 0 K is written as L 0 k .The index j denotes a conversion process between a source water species k s j and a target water species k t j .The effect of this process is expressed through the pseudoflux R j .We consider an arbitrary number m of such conversion processes.We now define variables λ kj = δ k,k s j − δ k,k t j for k = 1, . .., n and for j = 1, . .., m, where the usual definition of the Kronecker delta is used.The variable λ kj takes a value of 0 if a species k is not involved in the conversion process j ; it takes a value of −1 if it is the target species of this process; and it takes a value of 1 if it is the source species of this process.The variable λ kj will allow to write the time tendency of a water species by summing over all conversion processes, regardless of the role this specific water species plays in each process.Furthermore, a variable 0 This variable is the latent heat released at temperature T 0 under a conversion process with source species k s j and target process k t j .To clarify these notations, consider the original system of CGTBT07 with five water species and six conversion processes between them.By convention, we assign k = 1, . .., 5 to water vapour, liquid cloud water, precipitating rain, cloud ice crystals, and precipitating snow, respectively.Tables 1 and 2 give the values of λ kj and 0 j for the different conversion processes.
Next, a precipitation flux P k is defined for each component, even for the non-precipitating species (dry air, vapour, liquid cloud water droplets, and cloud ice crystals).Contradictory as this may sound, it should be stressed that in our barycentric system, these fluxes express the motion of the species with respect to the center of mass of the parcel.When Table 1.Variables λ kj for the system of CGTBT07.
precipitating species are present, the suspended species will move upward with respect to the mass center.Using a similar calculation as before to describe the motion with respect to the center of mass of the parcel, the relative precipitation fluxes P k are determined from the absolute fluxes P * k as where the absolute precipitation fluxes of suspended species can be taken to be zero.It should be noted that the strict distinction in CGTBT07 between suspended and precipitating species is somewhat arbitrary and scale dependent.Indeed, also the so-called suspended cloud water species can undergo a slow sedimentation.This arbitrary distinction is no longer necessary in the generalized set of equations that is presented here.Similarly to defining (relative) precipitation fluxes for all species, also diffusive fluxes D k are defined for all species, where the diffusive fluxes of precipitating species can be taken equal to zero.These notations make it possible to formulate the specific mass equations and the thermodynamic equation as follows: These equations generalize the ones from CGTBT07 in three ways: (i) an arbitrary number n of water species is considered; (ii) an arbitrary number m of interspecies conversion processes is considered; and (iii) the strict distinction between suspended and precipitating species can be abandoned.The fact that quite compact equations are obtained, which are valid for all components of the atmosphere, is an additional indication of the strength of the barycentric approach.

Remarks
Some comments should be given on the application area of the physics-dynamics interface presented in Eqs. ( 12)-( 13).
-The fact that these equations are very general, opens the road for a "plug-compatible" view of physics parameterizations.Indeed, the only output that is needed from a parameterization are diffusive and precipitative transport fluxes, pseudofluxes for phase changes, and the radiative and diffusive energy fluxes.The physics-dynamics interface then receives these quantities and determines the effect on the prognostic variables of the model, thereby ensuring satisfaction of the conservation of mass and energy, as well as consistency in the thermodynamic assumptions.However, it should be kept in mind that other conditions should be met before parameterizations can really be considered plug-compatible.A first aspect is that interactions exist between parameterizations.For instance, the parameterization of cloud processes will affect the radiation scheme.These kinds of interactions should properly be accounted for when plugging a new parameterization into a model.In this context, it is interesting to see that the technical recommendations that were made in Kalnay et al. (1989) regarding the design of parameterizations and their interactions, are still relevant at present.A second aspect is that parameterizations should also obey the second law of thermodynamics (Gassmann and Herzog, 2015).This condition cannot be enforced at the higher level of the physics- dynamics interface, and should be taken care of at the level of the parameterization itself.
-A common assumption in atmospheric modelling (although it is often made implicitly) is that all vertical mass transport due to the physics parameterizations is compensated for by a fictitious flux of dry air (Courtier et al., 1991).This assumption ensures the conservation of total mass in the atmosphere, but makes it impossible to express a net mass exchange with the surface due to, for instance, precipitation.From a barycentric point of view, this approximation means that the center of mass of an air parcel does not move vertically.The Eqs. ( 12)-( 13) remain valid under this assumption, if the absolute flux of dry air is defined as -The Eqs. ( 12)-( 13) are theoretically only valid for a model using the hydrostatic primitive equations.In a fully compressible system, the diabatic heating from the physics parameterizations does not only affect the temperature equation but also the continuity equation (Laprise, 1998).CGTBT07 present the extension of their flux-conservative system to the fully compressible case.An entirely equivalent development can be made for the generalized equations presented in this paper.
However, as shown by Malardel (2010), the impact of including the heat from parameterizations as a forcing in the continuity equation is quite limited.In other words, one can apply the thermodynamic Eq. ( 13) also in a nonhydrostatic model.
-The fact that Eq. ( 13) describes the evolution of enthalpy h = c p T , does not mean that this variable should become the prognostic thermodynamic variable of the model.A model that uses temperature T as the prognostic thermodynamic variable, can also use the presented interface.After all, one can easily calculate the total heat capacity tendency as follows: which in turn can be used to determine the temperature tendency from the enthalpy tendency: The importance of writing Eq. ( 13) as a time evolution of enthalpy only becomes clear in the time-discretized case.
where a superscript t denotes variables at the current time step, while a superscript t + t denotes variables at the next time step.Using an enthalpy-based formulation of the interface is reflected in the use of c t+ t p on the right-hand side of Eq. ( 16).Although this appears to be a small detail, it is crucial in ensuring the conservation of energy.The importance of appropriately discretizing a conserved nonlinear variable such as enthalpy is also indicated by Gassmann and Herzog (2008).
As a side remark, it can be noted that simply adding temperature tendencies from several parameterizations cannot lead to an energy-conserving atmospheric model, at least not for a process-split coupling strategy (Williamson, 2002).For example, consider a model containing two parameterizations (indicated with a and b), yielding a respective change in temperature of T a and T b , and a respective change in heat capacity of c a p and c b p .Suppose that each of these parameterizations is energy conservative in itself, meaning that the enthalpy changes are, respectively, h a = (c t p + c a p )(T t + T a ) − c t p T t and h b = (c t p + c b p )(T t + T b ) − c t p T t .Then the joint effect of the parameterizations cannot be expressed as T = T a + T b , but it should be determined as from which the total change in temperature is determined as This expression is only valid for a process-split coupling.For a time-split coupling, the total enthalpy change is still equal to the sum of the enthalpy changes of the separate processes, as indicated in Eq. ( 17).However, the fact should be taken into account that process b does not start from c t p and T t , but rather from the atmospheric state after accounting for process a, i.e.Working out the heat capacity and the temperature at the end of the time step now gives So with time-split coupling, the total temperature change can be obtained as the summation of the temperature changes from the separate parameterizations.However, it is better to use an enthalpy-based system, as this works both for the process-split and the time-split cases.
-The Eqs. ( 12)-( 13) only describe the evolution of the atmospheric prognostic variables.The prognostic variables of the surface scheme are not part of this system.In this context, the work of Best et al. (2004) should be mentioned.They present a method to separate the surface scheme from the atmospheric model.The core of this method is to describe the interaction between atmosphere and surface with fluxes.In this sense, their work matches perfectly with the flux-based Eqs. ( 12)-( 13).
3 Application of the flux-conservative equations in the AROME model AROME is a limited area model that was developed at Météo-France and is now a configuration inside the ALADIN system.It became operational in France in 2008, and it is currently used in many European countries of the ALADIN and HIRLAM consortia.AROME uses a nonhydrostatic, fully compressible dynamical core (Bubnová et al., 1995;Bénard et al., 2010), with the same spectral semi-implicit, semi-Lagrangian space-time discretization as the ECMWF's IFS model, and a terrain-following, mass-based vertical coordinate.(Laprise, 1992).AROME is coupled to the externalized surface scheme SURFEX (Masson et al., 2013) with the flux-based interface of Best et al. (2004).The physics parameterizations in AROME originate from the Meso-NH research model (Lafore et al., 1998).The Meso-NH model has a dynamical core which is explicit in time, with a staggered spatial grid and a height-based vertical coordinate, so it is substantially different from the AROME dynamical core.
The plugging of the physics from this model to a different dynamical core was quite challenging, and several approximations were made during this process.
A first approximation that is made in the existing AROME physics-dynamics interface concerns the heat transport by precipitation.From Eq. ( 13), it is clear that precipitation has two thermodynamic effects.Falling species modify the composition of the atmosphere, so they also change the specific heat capacity c p = n k=0 c k q k .Secondly, if a vertical temperature gradient exists, falling species are heated, thus cooling down the surrounding air.The effect on the enthalpy due to a change in c p is given by while the second effect due to a vertical temperature gradient is given by The combination of these two effects indeed corresponds to the effect of precipitation on the right-hand side of Eq. ( 13): The approximation made by the existing physicsdynamics interface in AROME is that it neglects the heat transport effect of precipitation, i.e. the term given in Eq. ( 19).
A second approximation concerns the effect of diffusive moisture transport (shallow convection and turbulence) in the energy budget.Similar to the effect of precipitation, diffusive moisture transport modifies the total specific heat capacity c p , and this effect should be accounted for in the energy budget.However, this effect is neglected in the existing AROME physics-dynamics interface.
www.geosci-model-dev.net/9/2129/2016/Geosci.Model Dev., 9, 2129-2142, 2016 2136 D. Degrauwe et al.: Generalization and application of the flux-conservative thermodynamic equations q q q q q q q q q q q 0 5 10 15 20 25 30 −0.5 0.0 0.5 1.0 1.5 Forecast range (h) MSL pressure (hPa) q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 2 m temperature (K) q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 25 30 Forecast range (h) 10 m wind speed (ms ) − 1 q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 25 30 −5 0 5 10 15 Forecast range (h) 2 m relative humidity (%) q q q q q q q q q q q Figure 2. RMSE (solid line) and bias (dashed line) over the period 1-30 November 2014, for REF (blue circles) and FCI (red triangles).
A third approximation is that the values of specific heat capacity c p and latent heat L i|l (T ) are not consistent between the different parameterizations.For instance, the heat capacity in the radiation scheme only accounts for water vapour and neglects the other hydrometeors (c rad p = (1 − q v )c pd + q v c pv ).This situation stems from the fact that the different physics parameterizations are developed by different teams, each using their own conventions.
A final approximation by the existing physics-dynamics interface in AROME is that the total temperature tendency is obtained by summing the temperature tendencies from the individual parameterizations.As indicated in the previous section, such an approach cannot lead to an energy-conserving system in a model with a process-split time step organization.
Although it can be expected that the overall effect of these approximations and inconsistencies is quite limited, the generalized physics-dynamics interface as presented in the previous section offers the possibility to get rid of them in order to take a (admittedly small) step towards a more accurate model.A second motivation to equip the AROME model with the generalized flux-conservative physics-dynamics interface is that this opens the route towards importing physics parameterizations from other NWP models, thus allowing a fair comparison of different parameterizations and stimulating scientific progress.

Impact on weather forecast
The impact of the presented flux-conservative formulation of the physics-dynamics interface is investigated with the AROME operational high-resolution LAM model running at Météo-France.Before April 2015, this model ran on a 739×709 grid with a resolution of 2.5 km. Figure 1 shows the model domain.The time step is 60 s.The model is provided with lateral boundary conditions by the operational global model "ARPEGE" from Météo-France.The initial conditions are generated with a 3DVAR data assimilation (Fischer et al., 2005;Brousseau et al., 2011).
At the surface level, precipitation and evapotranspiration imply a net mass flux across the surface.Since the vertical coordinate of the AROME model is mass based, correctly accounting for such net mass exchange between atmosphere and surface has far-reaching implications, especially in the surface boundary condition of the nonhydrostatic dynamical core.Currently, this has not been implemented in the dynamical core of the AROME model.Instead, the above-mentioned approximation is made that all vertical transport due to the parameterizations is compensated by a fictitious flux of dry air.Taking full advantage of the barycentric framework of Eqs. ( 12)-( 13) would require an adaptation of the dynamical core of AROME, which falls outside the scope of this work.
All these settings are identical for the operational run (denoted REF) with the temperature tendency-based interface and for the run with the flux-conservative interface (denoted FCI).

Monthly scores
The daily forecasts during two periods are considered in this section: 1-30 November 2014 and 6 January-6 Febru-Geosci.Model Dev., 9, 2129-2142, 2016 www.geosci-model-dev.net/9/2129/2016/q q q q q q q q q q q 0 5 10 15 20 25 30 −0.5 0.0 0.5 1.0 1.5 Forecast range (h) MSL pressure (hPa) q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 25 30 2 m temperature (K) q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 25 30 Forecast range (h) 10 m wind speed (ms ) − 1 q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 25 30 −5 0 5 10 15 Forecast range (h) 2 m relative humidity (%) q q q q q q q q q q q Figure 3. RMSE (solid line) and bias (dashed line) over the period 6 January-6 February 2015, for REF (blue circles) and FCI (red triangles).To avoid the problem of the double penalty, the precipitation is verified with the neighbourhood observation Brier skill score (Amodei and Stein, 2009).This score is determined by calculating the probability that a precipitation threshold is exceeded in the vicinity of an observation.By choosing the threshold, one focuses the verification more on light or on heavy precipitation.
The scores indicate that the impact of using the fluxconservative set of equations is quite limited when considering time-and space-averaged scores as the ones presented here.It should be stressed that no retuning has been done for the experiments with the flux-conservative equations.As a result, compensating errors can be responsible for masking an improvement of the scores.The fact that the scores do not change substantially, merely indicates that the approximations that are made in the existing temperature tendencybased interface are indeed small on a domain-wide scale.In this context, the limitations of this standard verification against station data should also be mentioned.By taking the average score over a large number of stations, important local differences may be hidden in the scores.In a similar way, the fact that monthly-averaged scores are considered, only allows to detect differences that are systematic in time.Therefore, notwithstanding the neutral impact on the standard scores, some significant differences may be observed under specific circumstances.A case study is presented in the next section to illustrate this.

Case study of a cold pool originating from heavy precipitation
When precipitation evaporates while falling through unsaturated air, it cools its environment.As such, a region of relatively cool air, the so-called cold pool, originates when heavy, localized precipitation occurs, for instance in precipitating convective systems (Fujita, 1959).It has been shown that the cold pool is in fact a key element in the life cycle of such systems.On the one hand, new convective cells originate at the border of the cold pool and its warmer surroundings, but on the other hand, if the cold pool becomes too strong, it may cut off the supply of warm air to the updraft (Engerer et al., 2008).The cold pool is also accompanied by a mesoscale high pressure area (Fujita, 1959) which plays a crucial role in the wind gusts that go with heavy precipitation.For these reasons, it is no surprise that an appropriate representation of the cold pool mechanism is essential in a NWP model (Engerer et al., 2008;De Meutter et al., 2014).Although evaporative cooling is the main cause for a cold pool, a second mechanism may enhance it.As precipitation falls from colder layers aloft to hotter layers below, it will be heated by the surrounding air, which in response will cool down (Johnson and Hamilton, 1988).As explained in Sect.3, this secondary thermodynamic effect (the transport of sensible heat) of precipitation is neglected in the existing AROME physics-dynamics interface, while it is correctly accounted for with the presented set of flux-conservative equations.One can thus expect that the intensity of a forecasted cold pool depends on which set of equations is used.This is confirmed when looking at the AROME forecasts over the Balearic islands on 19 January 2015.This case is characterized by convection developing ahead of an active cold front coming from the south.Figure 6a and b show the forecasted 12:00-18:00 UTC accumulated precipitation with the existing AROME interface (REF) and with the fluxconservative interface (FCI).It is observed that the overall structure of the precipitation is quite similar.However, when comparing the cold pool characteristics of both experiments, important differences appear.Figure 6c and d show the differences between both experiments for the 2 m temperature and the surface pressure.The temperature is significantly lower with FCI (up to 5 K cooler), and the surface pressure is higher (up to 1.4 hPa).
To further illustrate the impact of the heat transport by precipitation on the cold pool, the vertical profiles in the point as marked in Fig. 6b are studied for the experiment with the flux-conservative interface.The vertical profile of the precipitation fluxes (Fig. 7a) shows how snow and graupel originate aloft, they melt to form rain at around 850 hPa, and the rain starts to evaporate below 930 hPa. Figure 7b shows the vertical profile of the two phenomena that are responsible for the development of the cold pool, averaged between 12:00 and 18:00 UTC: the latent heat effects from phase changes (solid line), and the falling of cold hydrometeors into warmer air layers (dashed line).It is clear that the second effect is orders of magnitude smaller than the first effect, at least when considering the full vertical extent of the model.However, as shown in Fig. 7c, the heat transport by hydrometeors is not entirely negligible in the range between the surface and 900 hPa, and thus contributes to the intensity of the cold pool.
No comparison with observations is done for this case, because the purpose of this case study is merely to illustrate that even small terms in the energy budget can have a significant impact under certain conditions.The conclusions from this case study are in line with the results from Bryan and Fritsch (2002), where neglecting a supposedly small term in the energy budget unexpectedly leads to the worst results.

Conclusions
This paper starts from the equations presented in Catry et al. (2007) that describe how the effect of physical parameterizations on the dynamical core of an NWP model can be expressed in a flux-conservative way.The main advantage of these equations is that they impose the constraints of energy and mass conservation at a higher level in the model than at the level of the individual physical parameterizations.The presented equations only guarantee conservation of mass and energy regarding the effect of the physics contributions, not for the dynamical core of the model.A second advantage of the presented equations is that by gathering the thermodynamic calculations of all physics parameterizations in a single equation, it is also guaranteed that a predefined framework of hypotheses is consistently respected.
Notwithstanding these clear advantages, the equation set in the mentioned paper also faces limitations that hinder its application in existing NWP models.This paper presents a generalized set of thermodynamic equations that overcomes these restrictions without touching the sound theoretical foundations.More specifically, the presented equations are valid for an arbitrary number of hydrometeors, and can be applied in a model with an arbitrary number of conversion processes between these water species.This has allowed to use this set of equations in the AROME NWP model, which currently uses a physics-dynamics interface that makes some ad-hoc approximations.By moving to the generalized fluxconservative equations, the effect of these approximations can be studied.
Monthly verification scores show that the overall effect of introducing the flux-conservative equations in AROME is quite limited.There is no significant improvement or degradation of these scores.Given the mentioned theoretical benefits of the presented equations, this means that the presented work is a valuable advancement of the AROME model.Moreover, it appears that substantial differences may exist in specific cases.A detailed study of a heavy-precipitation case gives the example of the formation of a cold pool, which is an essential mechanism in the life cycle of a convective event.
As it appears, one mechanism that contributes to the formation of this cold pool is the heat transport by precipitation.This effect is neglected in the existing AROME physicsdynamics interface, while it is correctly accounted for in the presented flux-conservative set of equations.In this specific case, this leads to a different surface temperature and surface pressure within the cold pool.A more systematic study of the effect of heat transport on the life cycle of a cold pool is left for future research.In this paper, this case serves as an illustration of the importance of correctly accounting for supposedly small terms in the energy budget, something that is achieved with the presented set of thermodynamic equations.
Besides offering a direct improvement of the thermodynamic budget of the physics parameterizations of the AROME model, the presented set of equations also paves the way for interesting future research.Especially the impact of the heat from physics parameterizations on the continuity equation, and the effect of accounting for the net mass exchange between the atmosphere and the surface, are topics that deserve to be studied in detail.

Code availability
The used ALADIN codes, along with all related intellectual property rights, are owned by the members of the AL-ADIN consortium.Access to the ALADIN system, or elements thereof, can be granted upon request and for research purposes only.

D.
Degrauwe et al.:  Generalization and application of the flux-conservative thermodynamic equations -The precipitation fluxes P r and P s are relative to the (moving) center of mass of the parcel.They relate to the absolute precipitation fluxes P * r and P * s c p = c t p + c a p , and T = T t + T a .So for a time-split coupling, the enthalpy change of process b becomes h b = ( c p + c b p )( T + T b ) − c p T .

Figure 1 .
Figure 1.Operational AROME domain with a resolution of 2.5 km.The markers indicate the temperature stations used for the monthly scores.The dashed line indicates the area of the case study of Sect.4.2.

Figure 6 .
Figure 6.Case of heavy precipitation on 19 January 2015.The arrow and the marker in subfigure (b) indicate the location of the profiles of Fig. 7.

Figure 7 .
Figure 7. Vertical profiles at 18:00 UTC in the point indicated in Fig. 6b for the run with the flux-conservative interface.(a) precipitation fluxes: rain (black solid line), snow (red dashed line), and graupel (green dash-dotted line); (b) cold-pool-generating phenomena: latent heat effects due to phase changes (black solid line) and sensible heat advection (red dashed line); (c) same as (b) but focused on near-surface areas.