A model for moist convection in an ascending atmospheric column

This article presents a single‐column model of moist atmospheric convection. The problem is formulated in terms of conservation laws for mass, moist potential temperature and specific humidity of air parcels. A numerical adjustment algorithm is devised to model the convective adjustment of the column to a statically stable equilibrium state for a number of test cases. The algorithm is shown to converge to a weak solution with saturated and unsaturated parcels interleaved in the column as the vertical spatial grid size decreases. Such weak solutions would not be obtainable via discrete partial differential equation (PDE) methods, such as finite differences or finite volumes, from the governing Eulerian PDEs. An equivalent variational formulation of the problem is presented and numerical results show equivalence with those of the adjustment algorithm. Results are also presented for numerical simulations of an ascending atmospheric column as a series of steady states. The adjustment algorithm developed in this article is advantageous over similar algorithms because first it includes the latent heating of parcels due to the condensation of water vapour, and secondly it is computationally efficient making it implementable into current weather and climate models.


Introduction
The atmospheric models currently used in numerical weather and climate prediction consist of the classical compressible Navier-Stokes equations in a thin atmospheric shell around the Earth together with the laws of thermodynamics, and laws governing physical behaviours such as phase changes, radiation and surface fluxes. It is not feasible on today's (or the next decade's) computer architecture to computationally simulate this system at all scales, thus the governing equations are spatially and temporally averaged prior to being solved. This however means that weather and climate phenomena which occur on grid scales smaller and faster than that of the averaged system need to be incorporated in some 'average' sense via subgrid-scale parametrizations.
One such subgrid-scale phenomenon which needs to be incorporated in this fashion is moist cumulus convection. This is an important attribute of the weather because moist convection is responsible for large vertical fluxes of entropy, moisture, mass and momentum especially in the Tropics. These fluxes contribute to the tropical circulation which itself is a fundamental part of the climate system. Moist convection is also responsible for much of the severe weather in the extratropics. Currently convection can only be directly represented in numerical simulations with a horizontal grid resolution of at most 1-2 km (and often much less); note this is considerably smaller than the 10 km (or even more) spatial scale typically used in today's general circulation models (GCMs). In order to model the contribution of these small-scale processes on the larger scales, and to test new computational methods, it is useful to consider one-dimensional, single-column-atmosphere models, such as those developed in Yanai et al. (1973) and Tiedtke (1989), and references therein.
Numerical single-column models of the atmosphere are often used in climate prediction studies (Ghan et al., 2000;Sobel and Bretherton, 2000 and references therein) as they are simplified models in which the physics of processes such as radiation and convection can be explored. It is also well known by practising weather forecasters that the occurrence and intensity of convection can be relatively well predicted by studying the vertical potential temperature and moisture profiles together with knowledge of the external forcing (Met Office, 1993, chapter 4). These single-column models have formed a significant element of many theoretical moist convection studies (Holt, 1989;Shutts, 1994;Lock and Norbury, 2011). Lock and Norbury (2011) in particular performed an analytical study of the upward and downward motion of air parcels * in an incompressible single column, and presented a numerical algorithm using a Lagrangian-based hard convective adjustment step. The results of this algorithm were validated against analytical calculations. The numerical algorithm consisted of a discrete air parcel version of the Lagrangian problem, where the adjustment step acts to rearrange the parcels in a statically stable column (which is buoyantly stable and nowhere supersaturated). This rearrangement of parcels is designed to model the underlying physics of convection, without attempting to explicitly solve for the detailed three-dimensional turbulent dynamics in the column. Instead the rearrangement models the unresolved shorttime-scale dynamics by an instantaneous transport of air parcels.
In numerical weather and climate simulations the vertical profiles of the thermodynamic variables need to be adjusted at each time step to ensure the atmosphere is stable (or stable enough) for subsequent time-stepping. Currently an algorithm based upon the moist available potential energy (MAPE), sometimes referred to as the available potential energy (APE), is frequently used in the adjustment phase, with the atmosphere modelled as a single plume (Randall and Wang, 1992;Wang and Randall, 1994). However, satellite observations suggest that cloud layers often form at multiple levels, which implies that the unresolved circulations can be considerably more complex than a single plume. Another algorithm used is based upon the convective available potential energy (CAPE) which considers the available energy for a single infinitesimal parcel as it ascends the column (Emanuel, 1994;Riemann-Campe et al., 2009). The rearrangement algorithm proposed in this article allows saturated parcels both to condense water vapour and to find their level of neutral buoyancy in the column using physical conservation laws, and is therefore a first step towards modelling these more complex vertical circulations. Mathematically, moist convection is very tricky to study due to the nonlinearity of the moisture-heat exchange process and both the magnitude and nonlinear dependence of the moisture released by a parcel upon the potential temperature and the pressure level (or height). In this work we model this nonlinearity in a step-like fashion, while in reality this nonlinearity is a rapid, but smooth, process due to other physical features such as local mixing, condensation and evaporation of water occurring on finite, but small, space and time scales.
In this article we consider a similar approach to Lock and Norbury (2011), except we develop a more efficient algorithm which includes compressible effects, and also allows us to check that the convective inhibition criterion (CIC) is not violated. This criterion states that an ascending air parcel requires sufficient buoyancy to move above another parcel at the level at which they meet and switch places. Here we reformulate the Eulerian description of the single-column model and show that the convection of moisture in the column can be considered using conservation laws. In particular our model conserves the moist potential temperature θ M := θ + Lq of each parcel in the adjustment phase for a non-radiative system, where θ is the potential temperature, q is the specific humidity and L is a (large physical) constant. Note that θ M is a linearization of the realistic empirical moist potential temperature which is given by equation (3.8.4) in Gill (1982). Since θ M is conserved on parcels of atmosphere as they move in the column, any moisture condensation is matched by a latent heat release. A parcel will condense moisture and release latent heat if it becomes saturated and it is able to rise to a point of neutral buoyancy above this initial point. A parcel * In this article, we follow the usual understanding that a parcel is made up of a large but finite number of material particles. Each parcel has homogeneous properties and is transported by the velocity field. Note that moisture, i.e. the amount of water vapour in a parcel, is treated as a property rather than as actual material water molecules. becomes saturated/supersaturated if q ≥ Q sat (θ , p) where Q sat is the saturation specific humidity given by the Clausius-Clapeyron relation, p is a pressure-height variable and t is time.
Our proposed algorithm, when applied to such (unstable, saturated) situations, will instantaneously rearrange the parcels to a stable equilibrium configuration while obeying the relevant conservation laws for air mass, thermal energy and moisture mass. By simultaneously rearranging all the parcels we neglect the effect of individual parcels moving during the adjustment phase, i.e. we have two states, the initial pre-adjusted state and the final adjusted state. Such an instability of a moist midlatitude air column can be triggered, for example, by the column ascending due to an extratropical weather system penetrating under and lifting the column. In the latter part of this article we incorporate this ascending column effect into the model by including a timedependent lifting term into the saturation specific humidity for the column. A related algorithm to the one considered here is given in Cheng et al. (2017), but unlike here they consider incompressible columns and the effect of sorting individually ascending parcels, as opposed to the simultaneous sorting of all parcels carried out in this article. Cheng et al. (2017) provide a proof for the well-posed convergence of the solutions of their algorithm to weak solutions of the original partial differential equation (PDE) problem for the time-dependent ascending column in the limit as the time-step t → 0. A key condition for the convergence of their algorithm is satisfying the convective inhibition criterion (CIC). The CIC is incorporated into the algorithm presented in this article to aid its convergence in the time-dependent ascending column problem. The algorithm considered here is not identical to that in Cheng et al. (2017) (but the algorithms do coincide in the limit t → 0), and is designed to be more computationally efficient to make it implementable into current weather and climate models.
In order to simplify the analysis while still including the main physical feature of moist convection, that is the moisture release and the consequent latent heat absorption, we assume that the phase change is between water vapour and liquid water only, where any ice phase is neglected. We also exclude the role of the clouds, and mixing within the column, as well as the entrainment and detrainment of air into and out of the column.
The process of rearranging atmospheric parcels to find stable solutions is also used in variational problems posed for the atmosphere, such as determining the MAPE of a column (Lorenz, 1978;O'Gorman, 2010O'Gorman, , 2011Stansifer et al., 2017). In these problems an integral functional is minimised with respect to possible rearrangements of the air parcels. In this article we find, for the moist convection problem, an equivalent (but extended) variational formulation to that of the algorithms mentioned above, and show that numerical results for minimising the functional agree with those of our adjustment algorithm.
The article is laid out as follows. In section 2 we highlight the moist convection problem in the (usual) Eulerian PDE framework, before converting to the Lagrangian mathematical framework in section 3. In section 3 we also show how the moist convection problem for discrete saturated air parcels can be considered as the conservation of the moist potential temperature θ M , under a mass-conserving rearrangement of air parcels in the column. In section 4 we describe a discrete numerical algorithm to solve the moist convection problem. In particular, this algorithm determines the equilibrated state of a buoyantly unstable moist atmosphere. In this section we also formulate an equivalent problem as the minimisation of an integral functional. Results of the algorithm for more general initial data are presented in section 5 for both the adjustment of an unstable static atmospheric column, and an initially stable but ascending column which becomes unstable. Conclusions and discussion are given in section 6.

Eulerian mathematical formulation
We consider a single-column model of the atmosphere at a given longitude and latitude confined between fixed pressure levels p 0 Environment P seudo − adiabat p = p0 p = ptop Figure 1. Schematic of ascending air in the pseudo-adiabat within a region of descending air in the environment. In the actual column it is assumed that the pseudo-adiabats are uniformly distributed throughout the column. and p top . The column is assumed to be in hydrostatic balance, ∂p ∂z = −ρ(T, p)g, where for an ideal gas the density ρ(T, p) = p/RT where T is the absolute temperature, g is the gravitational constant, and the z-axis points upwards, perpendicular to the Earth's geoid surface. In practice this ideal gas law is formulated with a moisture correction term (Gill, 1982), but the version used here incorporates the leading-order compressible effects. In the z variable the column occupies the space z ∈ [0, z top ] where z top designates a rigid lid above the region where convection is occurring. In this article we work with pressure as the independent (vertical) coordinate to aid with the mass conservation of the proposed model. A detailed discussion of this is given in section 3. The column itself is considered to be made up of two types of connected regions, one consisting of upward motion (known as the pseudo-adiabat) and one consisting of downward motion (known as the environment), see Figure 1.
As the exact relative sizes of the pseudo-adiabat and environment regions are typically unknown, it is more convenient to work in a regime where the respective upward and downward motions in each region are considered simultaneously. This approach was considered in Lock and Norbury (2011) and thus forms the starting point for our analysis.
In pressure coordinates the conservation of thermal energy and moisture, with the horizontal velocities set to zero, are given by where θ and q represent the potential temperature and specific humidity of the air respectively, defined as The constant R = 287 J kg −1 K −1 is the gas constant for dry air, c p = 1004 J kg −1 K −1 is the specific heat at constant pressure for dry air and ρ v is the mass of water vapour per unit volume of moist air. The one-dimensional material derivative is defined by Here ω is the time derivative of the parcel position p = p A (t), as defined later in Eq. (7), where A is the Lagrangian label for the individual air parcels in the column. Then ω = −ρ(T, p)gw where w is the vertical velocity of the individual air parcels. Note that the thermodynamic Eqs (1) and (2) are just that, thermodynamic; they represent the vertical transfer of heat energy and moisture only, and do not model vertical momentum without a horizontal coordinate. Thus ω needs to be specified. Lock and Norbury (2011) shows both numerically and analytically that the Boussinesq version of Eqs (1) and (2) contain both the upward branch and downward branch of air motion in a one-dimensional column, and the reader is referred there for more details.
The quantities on the right-hand side of Eqs (1) and (2) are source terms for the potential temperature and humidity, and include an imposed net radiative heating or cooling term r(p), and Q L (to be defined later in Eq. (6)) the latent heating associated with the release of water due to the condensation of the water vapour (into cloud, rain, ice or a combination of these phases). The physical constant L = L v /c p = 2490 K, where L v is the latent heat of vaporization, the conversion factor required to find the increase in potential temperature of an air parcel if all its (warm) vapour is converted into condensed water at the same temperature and the associated heat energy released back to the air.
In the proposed model it is assumed that there are no fluxes of mass, moisture or thermal energy at the top or bottom of the column. As relative buoyancy of a particle is a monotonically decreasing function of p, then the column is said to be statically stable if throughout the column, where the equals sign denotes a neutrally stable column. The system of Eqs (1)-(4) is similar to the tropical models for bulk cloud microphysics, where the vertical momentum equation is replaced by the stability condition Eq. (4) on the longer time-scales for which these models apply (Biello and Majda, 2010). The independent variables in Eqs (1) and (2) are time t and pressure p. However for plotting purposes, it is convenient to convert our pressure variable to a height variable. Thus we note that an ideal gas in hydrostatic balance satisfies which upon integrating gives In Eqs (1) and (2) we model the latent heating process using a step-like nonlinearity hidden inside Q L , which takes the form where Q sat (θ , p) is a known function denoting the saturation specific humidity given by the Clausius-Clapeyron relation (Emanuel, 1994;Goldman, 2008;Lock and Norbury, 2011), which is defined in Eq. (16). Essentially Eq. (6) states that if a parcel of air is supersaturated/saturated (q ≥ Q sat ) then it has the potential to convert some of that moisture into latent heat, which could cause it to rise within the column, so while Eq. (4) is sufficient for stability, there could exist an alternative rearrangement of the parcels, found by lifting parcels and condensing moisture, which also satisfies Eq. (4) with q ≤ Q sat everywhere. In this model we assume that the condensed moisture is immediately precipitated and hence plays no further role in the model. The movement of potential temperature and specific humidity in the column could be computed by integrating the full Navier-Stokes equations (or a reduced form of these) in both the environment and pseudo-adiabat regions with entrainment and detrainment fluxes moving mass between the two regions, using discrete methods such as finite volumes etc. Instead however, we choose to replace the complicated turbulent processes of adjustment in a three-dimensional sense by the bulk vertical transportation of individual air parcels. In this adjustment process we do not attempt to model the dynamics in the column; instead we leave them unresolved and assume them to be an instantaneous response. In this case, parcels can move a finite distance in zero time, consistent with the time-scales observed in the atmosphere. The key differences between the algorithm developed in this article and mass-flux convection schemes, such as those of Tiedtke (1989) and Gregory and Rowntree (1990), are: the algorithm makes no attempt to solve for the dynamics of the ascending parcels, all parcels instantaneously jump to their new level and instability of the column is instantly removed; only saturated parcels can rise up the column; mixing with stable air is neglected; and finally there are no assumed entrainment rates, only the actual values of θ and q are used. Note, in Eq. (6) the physics is simplified to exclude the feedback from the clouds or rain that form, and so does not consider the ice phase of the water. However, the effect modelled by Eq. (6) captures the leading-order effect of moisture condensation in the atmosphere, namely the associated release of latent heat, and its consequences for parcel stability and movement.

Lagrangian mathematical description
Our goal is to find stable steady equilibrium solutions to Eqs (1)-(6) which can be achieved by considering the column as a one-dimensional array of individual air parcels and writing the system in Lagrangian form (Emanuel, 1994;Bokhove and Oliver, 2006). We split the column into equal air mass sized pressure parcels, because in our compressible model a simple rearrangement of these parcels in the adjustment process conserves the air mass of each parcel. The equivalent rearrangement process in z coordinates (such as that in Lock and Norbury (2011)) should allow the parcels to expand as they ascend the column and contract as they are forced downward. This process is automatically captured in pressure coordinates.
From Eq. (6) we have two distinct situations to consider when describing how air parcels move within the column: either a parcel is unsaturated (q < Q sat ) or it is supersaturated/saturated (q ≥ Q sat ). For the convection scenarios considered in this article we assume that the net radiative cooling occurs on a time-scale much slower than the movement of the air parcels so that this effect is neglected; we set r(p) ≡ 0 in Eq. (1) to simplify our presentation.

Unsaturated parcels
For parcels with q < Q sat , Eqs (1) and (2) become which can be integrated to where p = p A (t) denotes the position of the parcel with Lagrangian label A at time t. These equations show that, in this case, parcels simply move in the column with their initial potential temperature and moisture values, i.e. they move adiabatically as if they were dry air parcels. In our model these unsaturated parcels only move upward if Eq. (4) does not hold at some pressure level. During the latent heating phase these parcels can only be forced downwards by ascending moist parcels passing through them.

Saturated parcels
For all parcels, including saturated parcels with q = Q sat , the combination Eq. (1) + L Eq.
(2) leads to The quantity θ M A := θ A + Lq A denotes the moist potential temperature for a parcel A and is a conserved quantity as we move upwards with the moist parcel along the curve in the (p, θ )-plane such that θ + LQ sat (θ , p) = θ M A for each parcel A. We denote this curve C(A). Thus any loss of moisture by a parcel must be accompanied by an increase in potential temperature through latent heating. In this work we make the assumption that all condensed moisture is removed from the column immediately (as rain), and consequently that descending parcels conserve q and θ , i.e. they do not re-absorb moisture.
The saturation specific humidity function Q sat is a smoothly differentiable function of (θ , p) which satisfies the conditions for all θ and p ∈ [p 0 , p top ]. See Eq. (16) for an explicit form for Q sat which is relevant to most of the Earth's atmosphere. Any saturated parcel with Lagrangian label A will remain on C(A) (i.e. it will lose water vapour, but will remain saturated) as it moves upwards. Along this curve (often referred to as the moist adiabat) the potential temperature θ A of air parcel A, with moist potential temperature θ M A as it rises to level p A , is given by the solution of the implicit equation

Moist convection adjustment algorithm
For a column of parcels the system is buoyantly unstable if Eq. (4) is violated anywhere in the column. If Eq. (4) does hold, the atmosphere is buoyantly stable, but if some of those parcels are saturated there could exist an alternative arrangement of the air parcels to which the column could adjust, in which Eq. (4) again holds. In that case the saturated parcels could rise to a new buoyantly stable level. In either of these situations the air parcels are convectively adjusted into a stable column by rearranging the positions of the parcels (via the unresolved dynamics). In the case of a dry column (q = 0 or q < Q sat (θ , p) for all parcels for all times), the adjustment of the air column is unique and θ is arranged monotonically increasing from the bottom of the column (Shutts and Cullen, 1987). However, for the moist column there are potentially many stable rearrangements of the parcels which satisfy Eq. (4) and have q ≤ Q sat everywhere.
Here we formulate a moist convection adjustment algorithm to adjust the column of air parcels subject to the conservation laws Eq. (7) or Eq. (8) depending on whether q < Q sat or q = Q sat respectively in each parcel. This algorithm has similarities to the top-down algorithm discussed in Wong et al. (2016) and Saenz et al. (2015) and the greedy algorithm presented in Stansifer et al. (2017). The novelty of our algorithm is that we propose to model the moist convective adjustment in a single column using a parcel rearrangement which takes into account the latent heating (from the condensing moisture) of each parcel as it rises. We then apply the convective inhibition criterion (CIC) which states that an ascending parcel requires enough buoyancy to move past another parcel at the level at which parcels meet and switch places.
The moist convective adjustment algorithm provides a parametrization of the latent heating processes in the mathematical model of section 3 where adjustment acts to ensure that the column atmosphere is statically stable (Eq. (4)), while parcels remain nowhere supersaturated (q ≤ Q sat ). We start by dividing the column atmosphere into N equal discrete pressure parcels withp for i = 1, . . . , N, so each parcel has equal mass p/g where p = (p 0 − p top )/N. Note that here the subscript denotes an Eulerian position of each parcel not the Lagrangian label of the parcel; we use the tilde over each variable to denote it is Eulerian. Thus at some time t, the pre-adjusted (or initial) state of the parcels at levelsp 1 ,p 2 , . . . ,p N have potential temperaturesθ 0 1 ,θ 0 2 , . . . , θ 0 N , moisture contentsq 0 1 ,q 0 2 , . . . ,q 0 N and Lagrangian position labels A 0 1 , A 0 2 , . . . , A 0 N . In the adjustment context, the Lagrangian position labels are used to track which parcels have moved to which new location in the column during the adjustment process. The value of the potential temperature and the specific humidity are assumed to be evenly distributed throughout each parcel and the parcels pass through each other during the adjustment without modifying either of these quantities. If our initial configuration is one where Eq. (4) does not hold, we first rearrange the column such that θ (p) is monotonically decreasing in p, i.e. such that Eq. (4) holds but with some parcels now supersaturated (q > Q sat ). We then apply the algorithm below.
-START -In the moist convection adjustment algorithm the superscripts 0 denote the initial pre-adjusted states, the superscripts 1 denote the final adjusted state, the superscripts U denote initially unsaturated parcels, and the superscripts T denote a temporary value while the sort is occurring. The moist convection adjustment algorithm proceeds from the top of the column as follows. For the 1. Temporarily move all unsorted parcels from their current level to this new sorting level for j = {1, . . . , N}\ {sorted parcels}. 2. (a) If j < N sort andq 0 j ≥ Q sat (θ 0 j ,p j ) (i.e. the parcel is saturated or supersaturated) then excess moisture will condense and its potential temperature will increase by the latent heat exchange. Thus its new potential temperature θ j is found by solving the implicit relation The hat on θ j denotes its Eulerian value during the adjustment sorting phase. The corresponding moisture level is then found via the parcel is unsaturated) then the parcel should not move. (However, for the sake of coding the algorithm one can assign θ j =θ 0 j asθ 0 j <θ 0 N sort so the parcel currently at level j will never be installed at level N sort .) (c) If j ≥ N sort andq 0 j ≤ Q sat (θ 0 j ,p j ) then these parcels move down adiabatically so θ j =θ 0 j and q j =q 0 j . (d) If j ≥ N sort andq 0 j > Q sat (θ 0 j ,p j ) then the excess moisture is converted into latent heat at level j and then the parcel is moved adiabatically to level N sort . Thus θ j is a solution to and q j is given by the first equation in Eq. (11).
3. The process in step (2) generates an array of potential temperatures (N sort ) = { θ k : k is an unsorted parcel} temporarily moved to level N sort . The largest of these potential temperatures, θ J fromp J say, is then checked to see whether it satisfies the CIC. If this parcel is unsaturated we are done.
If it is saturated we check that its potential temperature is greater than that of all the unsaturated parcels it had to pass to get to level N sort . Mathematically this means checking that where θ T J is a solution of for allp U k k ∈ (J, N sort ). If Eq. (12) holds for all k ∈ (J, N sort ) then the parcel from pressure levelp J is installed at the levelp N sort with The parcel from levelp J is then eliminated from the sort, and the algorithm moves to levelp N sort −1 . If Eq. (12) does not hold for all k in this range then this parcel cannot rise top N sort and we then check the CIC for the parcel with the second largest potential temperature in (N sort ) at the current sorting level. This process continues until either the parcel is unsaturated or Eq. (12) holds for all k ∈ (J, N sort ).
4. Steps 1-3 are repeated at each decreasing level N sort < N until N sort = 1 and every parcel is assigned a level. Then the sorted (adjusted) values areθ 1 i ,q 1 i , A 1 i for i = 1, . . . , N.
-END -Because all saturated parcels are temporarily ascended to every height until they are sorted (with the algorithm choosing the largest θ value satisfying the CIC for the post-rearrangement θ at that level) the algorithm produces a rearrangement which generates as much latent heating as the initial conditions will allow, thus condensing as much moisture as allowed. Numerical simulations of this algorithm for various test cases are presented in section 5. Solutions of the algorithm are shown to converge to either smooth solutions with jumps (here we mean smooth solutions in adjoining regions), or non-smooth solutions where saturated and unsaturated parcels are interleaved.
In the case of time-dependent convection problems (such as those considered in section 5.2) where the above algorithm is applied at each time-step, a mathematical proof for the convergence of a related algorithm to that presented here, as the time-step t → 0, can be found in Cheng et al. (2017). They found that the CIC was essential for proving well-posedness, and to keep the rate of mass transfer bounded as t → 0. As the above algorithm sorts all the parcels simultaneously it is difficult here to check this criterion exactly as the point at which a rising parcel meets other parcels is unknown. However, we make the assumption that if a saturated parcel needs to 'overtake' another saturated parcel as they both rise then they meet at a height where this overtaking can happen. Therefore we only check that the saturated parcel can move past all the descending unsaturated parcels; see step (3) in the numerical algorithm. As we sort all the parcels simultaneously we make the assumption that the saturated parcels meet the unsaturated ones at their initial preadjusted positions. Note that although the algorithm considered here is not identical to that in Cheng et al. (2017), as t → 0 the two algorithms coincide. This is because in this limit parcels will ascend one by one as they hit their respective descending C(A) curves.

Variational formulation
The adjustment algorithm developed in section 4.1 can be considered in an equivalent way as the minimisation of an integral functional. The functional in question takes the form where a > 0 is a constant. Here K(p, p ) is a so-called transport plan, defined on [p 0 , p top ] 2 which is unknown and needs to be determined such that it satisfies the relations and θ * (p, p ) is a function which gives the potential temperature of a parcel at level p, which prior to the rearrangement was at p . The transport plan K(p, p ) can be interpreted as the probability for a parcel at p to jump to p. This transport plan formulation has been used extensively in atmospheric science problems (Romps and Kuang, 2011;Stull, 2012;Stansifer et al., 2017), where it is also known as a transilient matrix. The two integral conditions on K(p, p ) in Eq. (14) ensure that air mass is conserved in the rearrangement, and the exponential in Eq. (13) weights the value of Kθ * such that more emphasis is put on the top of the column. In section 5.1 we perform numerical minimisations of F a [K], hence we consider Eq. (13) for a discrete set of fluid parcels. Our numerical minimisations of Eq. (13) show that for the discrete problem, K(p, p ) is a permutation matrix with only one non-zero entry (a one) in each row and each column. Our minimisation reorders the parcels, it does not mix them, which would correspond to K having entries which are not 0s or 1s. This fact is non-trivial but follows from Birkhoff's theorem (Birkhoff, 1946). This suggests that in the continuous limit (i.e. the number of discrete parcels in the discretization of Eq. (13) tends to infinity), that where g(p ) is a map and δ(·) is the Dirac delta function; but in general the limit must be taken in a probability measure sense, and the 'weak' limit is a transport plan, not a map. The function θ * (p, p ) in Eq. (13) takes the form where (p, θ M (p )) is the unique solution (given p, p ) to The first row of Eq. (15) converts moisture to latent heat for a parcel which ascends the column; the second row deals with any supersaturated parcels which are forced down by ascending parcels, in this case any excess moisture is rained off and then the parcel is moved down adiabatically. The final row of Eq. (15) deals with unsaturated parcels or those not satisfying the convective inhibition criteria (CIC). The construction of F a [K] is such that in the limit a → ∞ the optimal solution minimising θ (p), on a parcel level by parcel level basis, will be chosen prioritizing the top of the domain and working downwards. We claim, a claim which is verified numerically in section 5, that minimising Eq. (13) with a → ∞ gives the same solution as the algorithm set out in section 4. Thus the minimisation process in the limit a → ∞ will yield a monotonic θ (p) from Eq. (13) and produce the global minimum of Eq. (13), subject to Eqs (14) and (15). Note that Eq.
(15) can be modified to include different physical processes.

Numerical test cases
In section 5 we present results for some test cases in order to demonstrate the adjustment algorithm of section 4.1. To present numerical results we require an explicit form for the specific humidity function Q sat . The form of Q sat is given by an empirical formula set out in Gill (1982) and the Appendix of Lock and Norbury (2011), where where Q 0 is a positive constant and e sat (T) is the saturation vapour pressure of water. To approximate the function e sat (T), we use the formula for the saturation vapour pressure of pure water over a plane water surface (Gill, 1982, p. 606) log 10 (e sat ) = where Q 1 , Q 2 , Q 3 and T 0 are positive constants. The values of Q 0 , Q 1 , Q 2 , Q 3 and T 0 are given in Table 1. For 233 ≤ T ≤ 313 K the values of e sat using the above formula agree with those in the Smithsonian Meteorological Tables to within 0.2%. Thus combining the two above formulae gives The form of Q sat in Eq. (16) will be used in section 5.1 to test the adjustment algorithm in moving from a prior arrangement of air parcels to an adjusted profile. Once we have investigated the adjustment phase we will go on to demonstrate the algorithm's potential for investigating time-dependent problems. The particular time-dependent problem we investigate is that of an ascending air column being lifted at a speed α as a consequence of a cold front of extratropical air penetrating under the column (Figure 2), effectively thickening the thin planetary boundary layer (PBL) which is assumed to exist at the bottom of the column. In this layer the flow field is strongly influenced by interaction with the surface of the Earth (Holton and Hakim, 2012).
We assume that the lifting of the column occurs on a slower time-scale than the moisture condensation and the adjustment of the column, but on a faster time-scale than the net radiative heating or cooling, which we earlier neglected. Therefore we can adjust the column to its new arrangement before lifting the column. The lifting of the column essentially amounts to On the left is a schematic of the air column in consideration above the planetary boundary layer (PBL). On the right is the same column which has been lifted by a cold air front penetrating below the column. making Q sat time-dependent by lowering the Q sat function, and re-adjusting the column, hence here we modify Eq. (16) to where P(t) is determined such that the bottom of the column is lifted at a speed α; see section 5.2 for more information.
In the adjustment stage we do not model the dynamics within the column, which are assumed to be an instantaneous response, and the parcels can move a finite distance in zero time. This is consistent with the time-scales observed in the atmosphere where air parcels in an active moist atmosphere in the extratropics can move from the surface to the tropopause in the space of approximately 60 min.

Results
The results presented in this article are divided into two scenarios. In section 5.1 we consider the adjustment of a static column in order to present results for the moist convection adjustment algorithm outlined in section 4.1 and demonstrate its numerical convergence as the number of parcels N → ∞. We also compare the results of the adjustment algorithm with solutions to the direct minimisation of F a [K] from section 4.2. These minimisations are numerically obtained by dividing the column into N equalsized pressure parcels, see Eq. (10), and using the Munkres (or Hungarian) algorithm (chapter 4 of Burkard et al. (2009)), as this algorithm always finds the global minimum of the functional. When minimising F a we find a = 0.007 to be large enough to give agreement with the adjustment algorithm. While this seems like a small value for a, it is in fact large when the column pressure is non-dimensionalised to p = 1 at the planetary boundary layer, where it corresponds to the value 700, hence the results of the Munkres algorithm minimising F a presented in this section are in the a → ∞ limit (which we refer to as F ∞ ).
Secondly in section 5.2 we consider an ascending column, as in Figure 2, as a series of steady states and examine how the moisture and potential temperature profiles vary with time. In all the results presented in this section we set the pressure p top = 11 250 Pa and p 0 = 10 5 Pa.

Dry convection
In order to test the algorithm in a familiar scenario we first consider the case of dry convection. For the initial condition we take the unstable potential temperature distribution and q 0 = 0. The initial distribution for θ 0 is given in Figure 3 as a function of z, where z(p) is given by (5). The initial potential temperature distribution does not satisfy Eq. (4) hence the algorithm of section 4.1 should adjust the column such that Eq. (4) is satisfied. The results of the algorithm are presented in Figure 4 for N = 10 000. In the case of dry convection we know that the adjustment for θ is unique and that the parcels are simply arranged in monotonically increasing θ order. The results in Figure 4 show that the solution from the adjustment algorithm rearranges to this unique solution, and the squares in Figure 4(a) show the solution minimising F ∞ [K], which agrees with the adjustment algorithm.

Moist convection
For the moist convection problem we first consider the unstable initial potential temperature distribution given by For the initial moisture content, we want to create an air column containing regions of saturated and unsaturated air, hence we set where The initial distributions θ 0 and q 0 plotted as a function of z are given in Figure 5.  In Figure 5(b) those saturated parcels at the base of the column contain the most moisture and hence have the greatest latent heat pool. Therefore during the adjustment process it is these parcels which are expected to ascend furthest up the column. Figure 6 plots the potential temperature θ and the moisture q distributions as a function of z after the moist convection adjustment. The pairs of panels display the results for N = 100 and 10 000 from top to bottom. The results show that the saturated parcels at the bottom of the column condense a proportion of their moisture and ascend the column to new heights throughout the column in the range 3500 m z 9000 m. We can see that it is the lower parcels which have ascended the column by considering Figure 7 which shows the final Eulerian position of each parcel compared to its initial position in the column. These figures show that it is parcels with label initially less than 11 and 1126 respectively which ascend in the adjustment phase and these are parcels with z 830 m prior to the adjustment.
The results in Figure 6 show that as these lower parcels rise they displace parcels down into a region of the column where they are no longer saturated. This gives a region at the base of the column containing unsaturated parcels, above which occurs a region of the column where the moisture distribution is nonsmooth. Note, there is also a small region at the bottom of the column where q is non-smooth and these are parcels which contained a large amount of moisture prior to adjustment but  Figure 5. Plot of (a) the initial potential temperature distribution θ 0 (z) and (b) the initial moisture distribution q 0 (z) given by Eqs (19) and (20) respectively. In panel (b) the dashed line gives the value of Q sat . The adjustment of this data leads to the weak solution in Figure 6.
were forced down the column by monotonically reordering θ before applying the adjustment algorithm. In the higher of the two non-smooth regions, ascending saturated parcels are interleaved with unsaturated parcels, producing the discontinuous jumping seen in Figures 6(b) and (d). As N increases, the moisture profile converges to a profile with saturated parcels interspersed with unsaturated parcels throughout the main part of the column. This 'weak' type of solution has also been observed in the oceanographic work by Hieronymus and Nycander (2015), who globally minimize an enthalpy function to find the Available Potential Energy. This type of weak solution would not be attainable by regular discrete PDE methods. Note, despite the non-smooth moisture profile, the corresponding θ distribution appears smooth. The corresponding tephigrams for the initial data in Figure 5 and the result after the adjustment phase from Figure 6 are given in Figure 8. Figure 6(c) also contains results from minimising F ∞ [K] directly, denoted by the squares. This minimisation is carried out with N = 2000 parcels but we only plot every 20 for clarity. This shows that numerically the solution to minimizing F ∞ [K] is in excellent agreement with the adjustment algorithm result. We also get excellent agreement with the moisture profile too when minimising F ∞ , but this result is not shown. Note that if a is not large enough in Eq. (13) then the solution generated by minimising F a [K] is not a stable solution and there will be regions of the column where dθ dp > 0. The major advantage of using the adjustment algorithm over the Munkres algorithm is that numerically it generates the solution in O(N 2 ) operations as opposed to the O(N 3 ) operations in the Munkres algorithm. This speed-up time allows for time-dependent ascending column calculations, such as those in section 5.2, to be carried out in time frames suitable for numerical weather prediction models.
The results in Figures 6 and 7 show that the adjustment of this initial condition is akin to the process which occurs in deep convective systems such as thunderstorms. In these systems the fluid parcels undergo a large scale adjustment allowing parcels at the base of the column to condense a large volume of moisture as they ascend to a new level of neutral buoyancy. This large volume of condensed moisture produces heavy precipitation events. In shallow convection systems this may not be the case and saturated parcels near the bottom of the column may be halted in their ascent at lower levels by reaching some alternative stable state. In the case of shallow convective systems the pre-adjustment profile will likely be different to that above, with a θ distribution more stable than that considered in Figure 5.
Despite the solution produced by the adjustment algorithm in Figure 6 being 'weak', we observe convergence as N → ∞. We show this in the numerical results of Figure 9(a). This figure displays the total mass of moisture content of the column after the adjustment as a function of N. Here the integral is calculated via the trapezoidal rule. Figure 9(b) plots the corresponding value of the initial mass of moisture prior to the adjustment phase q Tot 0 where the exact value is known (dashed line), to show the convergence of the trapezoidal integration scheme. This figure shows that for N 10 000 the solution has effectively converged, with the error in the value of q Tot 0 < 0.2% of the exact value at N = 10 000. As the computational run time of the moist convection adjustment algorithm scales ∝ N 2 , in order to efficiently, but accurately, calculate results for the time-dependent ascending column problem we need a trade-off between run time and accuracy. As we are interested in this article in identifying the nature of the weak solution we choose N such that the results are converged to graphical accuracy; N = 10 000 is a good choice for use in the time-dependent results in section 5.2. However, in global numerical weather prediction models it may be that a smaller value of N is sufficient and practical. The jumps in q Tot observed in Figure 9(a) are due to the weak form of the solution for q.
Not all solutions to the moist convective adjustment algorithm have rapid variation in the final moisture distribution, and in Figures 10 and 11 we show an example which leads to a nearly smooth solution post-adjustment. The initial data for this example has and This example has a moist layer of parcels at the base of the column, all of which rise to the very top of the column in the adjustment phase. This leads to a substantially hotter region at the top of the column which contains a small amount of un-condensed moisture. We term this solution a nearly smooth solution because θ and q are each smooth in two distinct adjoining regions z ≶ 11 500 m. Discrete PDE methods mostly produce a good approximation to this solution with some smoothing around the discontinuity; this is a classically weak PDE solution.
After this brief analysis of the moist convection adjustment algorithm, we use it to analyse the problem of an ascending vertical column of air.

Ascending column
The moist convection adjustment algorithm analysed in section 5.1 would typically be used at each time step in a numerical weather or climate model to ensure the atmosphere is everywhere stable (or stable enough) in order to perform the subsequent time step. Currently simulations at the UK Met Office undertake this stabilization using a convective available potential energy (CAPE) adjustment process. CAPE considers the available potential energy for a single infinitesimal parcel, ignoring any effect it might have on the environment region. In this section we perform a similar calculation by modelling the ascent of a column of air which is lifted at a speed α by a cold front of extratropical air penetrating under the column (Figure 2). This effectively amounts to lowering C(A) for each parcel at a speed α, which is achieved by the following choice of P(t) in Eq. (17). The function P(t) is defined (numerically) at the nth time step, t n , as P(t n ) = P(t n−1 ) P, where P(t 0 ) = 1, and P is the (interpolated) solution to The function z n−1 (p), from Eq. (5), is z(p) from the (n − 1)th time step. The above interpolation determines the pressure Pp 0 at the base of the column when it is lifted an amount α z, hence P(t n )p 0 is the result of n time steps and hence n lifts. This ensures that the base of the column rises at approximately the constant speed α, by lowering the Q sat function at the same rate. As the Q sat function is lowered, if, and when, some parcels become saturated or supersaturated, we then apply the adjustment algorithm. Note that we choose to lower each C(A) curve in this study; equally we could have modified our code to lift the air column. The results are identical, but this approach means the discretization of the column Eq. (10) is fixed.
For the results in this section we consider an air column ascending at the approximate speed α = 125 3 m h −1 and we take time steps t = 1 hour over the period of 96 h, with N = 10 000. As we adjust the parcels simultaneously, this algorithm effectively has a minimum time step t adj , which is the time-scale for a parcel to ascend from the Earth's surface to the bottom of the tropopause in reality. Therefore the assumption of instantaneous  Figure 10. Plots of (a) the initial potential temperature distribution θ 0 (z) and (b) the initial moisture distribution q 0 (z) given by Eqs (21) and (22). In panel (b) the dashed line gives the value of Q sat . The adjustment of this data leads to the nearly smooth solution in Figure 11. adjustment of the column will be a good approximation provided t adj t, hence the choice of t = 1 h. However, the algorithm works fine with t < t adj . We consider an initially stable column with and Thus the parcels are 90% saturated at p = p * , have a constant moisture content for p ≥ p * , and tend to 80% saturated at the top of the column. See Figure 12 for plots of Eqs (23) and (24). Note, we use the values of p * = p 0 , 89 150, 79 300, 70 380 Pa which give a corresponding value of z * = z(p * ) = 0, 1000, 2000, 3000 m. The time evolution of the total moisture level q Tot (t) is plotted in Figure 13 for these values of z * .
The results for q Tot (t) show that the z * = 0 m, z * = 1000 m and z * = 2000 m initial data undergo a large adjustment after a few hours when the column first becomes saturated at some  Figure 11. Plots of (a) the final potential temperature distribution θ(z) and (b) the final moisture distribution q(z) after the moist convective adjustment for the initial data given in Figure 10. The dashed line in panel (b) gives the value of Q sat after the rearrangement.
height. This large adjustment leads to a rapid condensation of moisture accompanied by an associated potential temperature increase in the column. After this large adjustment the decrease in moisture is more gradual, i.e. a more constant precipitation event. The z * = 0 m and z * = 1000 m data contain enough moisture to have additional, smaller, precipitation events between 12 and 36 h. These significant precipitation events are due to the sudden release of MAPE from the column, which was previously unavailable due to a lack of saturation. This displays the severe precipitation which convection is responsible for in some weather events. The z * = 3000 m data on the other hand does not undergo a large adjustment during the lifting of the column as the parcels at the base of the column do not contain enough potential latent heat to allow them to ascend up the column very far. In this case the column continues to condense moisture at an approximately constant rate over the ascent of the column. The size and positions of the large and smaller adjustment events are somewhat dependent on the values of N and t chosen, but in the limit t → 0 with N → ∞ we do see numerical convergence, where only one parcel ascends in each time step. However, for practical numerical weather prediction models, where N and t will not be in this limit, some grid dependence is to be expected.
The moisture plots in Figure 14 show that by t = 72 h all the profiles are almost completely saturated with a large increase in potential temperature at the base of the column. Beyond this time  Figure 12. Plot of (a) θ 0 (z) and (b) q 0 (z) for the ascending atmospheric column given by Eqs (23)  none of the parcels at the base of the column are able to rise up the column in a rapid adjustment process. In this period the moisture loss shown in Figure 13 occurs at a similar rate.
This is a simple model of an ascending column with moisture condensation and heat transfer; many additional features such as entrainment and detrainment of moisture and the radiative cooling of the column are neglected, but the results in Figure 14 show that the adjustment of a column atmosphere can lead to a large amount of condensation of moisture and hence a rapid increase in the potential temperature of the column. The novelty of our algorithm is that this large moisture condensation is directly related to the latent heat increase of an air parcel via the conservation of the moist potential energy. This potential temperature increase and moisture reduction will have a large knock-on effect on the dynamics of numerical weather or climate models, showing that accurately modelling this adjustment phase may be vital to these systems.

Conclusion and discussion
In this article we develop a mathematical approach for investigating moist convection in a compressible single-column atmosphere. The model captures the leading-order effects of moist convection (namely the conversion of water vapour to moisture which releases significant amounts of latent heat). The single column has no horizontal variation (or degrees of freedom), and determines steady, stable solutions of the dynamic atmosphere. The model is based upon conservation laws for the thermal energy in terms of the moist potential temperature θ M = θ + Lq, the air mass and the moisture mass, which are built into the governing Eulerian equations, and are more explicit in the Lagrangian approach.
The conservation laws are used to construct a numerical algorithm to model the moist convective adjustment of an unstable column of atmosphere. This algorithm is novel in that it allows saturated compressible parcels to find their level of neutral buoyancy using physical conservation laws, and it incorporates the convective inhibition criterion, which only allows parcels to pass through each other if the lower parcel has sufficient buoyancy. In this algorithm the adjustment of the column is assumed to occur instantaneously, hence the parcels can move a finite distance in zero time, and this is consistent with the time-scales observed in the atmosphere where air parcels can rise from the surface to the tropopause in approximately 60 min in an active moist atmosphere in the extratropics. So essentially there is a lower bound to the size of our algorithm time step in time-dependent adjustment problems such that the assumption of instantaneous adjustment is a good approximation.
The moist convective adjustment algorithm presented here adjusts the atmosphere to produce a rearrangement of the moist air parcels which minimises the functional F a [K] given in Eq. (13) by allowing each parcel to gain a (large) amount of latent heat from condensing moisture, subject to the system constraints. The numerical results of this adjustment algorithm for a moist column show, for appropriate initial conditions, saturated parcels interspersed with unsaturated parcels in the vertical direction, which we term a 'weak' solution, where this interleaving is preserved under mesh refinement. These 'weak' solutions are well-posed solutions to the moist convection problem (Cheng et al., 2017), but they would not be obtainable from the Eulerian PDEs using conventional discrete PDE methods such as finite differences or finite volumes, and hence different or new approaches are needed to deal with these types of 'weak' solutions in computational weather and climate systems. Despite the 'weak' nature of the solutions it is shown that the algorithm converges numerically as the vertical mesh size is refined. This algorithm could be used at each time step in a numerical weather or climate model to stabilize the atmosphere, as an alternative to current algorithms based upon the moist available potential energy (MAPE) or the convective available potential energy (CAPE).
The adjustment algorithm was also used to calculate results for an ascending column of atmosphere to mimic the case of a moist weather system vertically displaced by a cold extratropical weather front penetrating underneath it. In this problem parcels of atmosphere became saturated as the column was lifted, causing the atmospheric column to become unstable. The moist convection adjustment algorithm was then applied at each time step to adjust the column parcels. Here it was found that if the column contained sufficient moisture then the adjustment could cause a large amount of moisture to condense, implying large amounts of precipitation. When the column initially did not contain sufficient moisture then the column continues to condense moisture at a slower, more constant, rate during its ascent.
The adjustment algorithm presented here can be extended to include more physical effects (and hence to more accurately simulate reality); these effects include radiative cooling, mixing of temperature and moisture as parcels exchange positions, freezing of condensed water, and surface heating. This column model can be extended by including entrainment/detrainment of air between horizontally adjoined columns to simulate two-dimensional (2D) and 3D scenarios, as in Wong et al. (2016); this would add moist convection processes to their work. However, the inclusion of some of these additional features is non-trivial as additional timescales are introduced, and the constraint of a single column will limit the modelled reality in any case.