Modification of the thermodynamic variability closure in the Met Office Unified Model prognostic cloud scheme

In the Met Office Unified Model prognostic cloud scheme (PC2), thermodynamic variability is represented by considering two distributions, based on the properties of the clear and cloudy portions of the grid‐box. To calculate a cloud fraction increment, the equation set needs to be closed. The closure assumption determines how these two distributions are blended together at the saturation boundary. The current closure is reviewed. In the limit of vanishing liquid condensate, this can lead to very tall but narrow probability distribution functions (PDFs). At present, the two parts of the PDF are weighted based on the area under their curves. As a result, these very tall PDFs are allowed to contribute significantly, since their area is finite. This closure is shown to sometimes lead to unphysically large cloud fraction changes. However, they are prevented from affecting the model evolution by being overwritten if their magnitude is unrealistically large. A new blending closure is proposed. This weighs the two PDFs based on their widths. The contribution of very tall, but narrow, PDFs is hence drastically reduced. While this removes the undesirable need for overwriting unrealistic increments, it is shown that overall model performance did not improve or worsen.


| INTRODUCTION
While clouds play an important role in the Earth's hydrological cycle and radiation budget, the processes that determine their formation and dissipation occur on scales much smaller than those that can be resolved by the global numerical models used for performing weather forecasts or climate projections on the global scale. As a result, atmospheric general circulation models make use of cloud cover parameterisations (Stensrud, 2007). These are schemes that calculate the cloud cover and the condensate amount in each model grid-box at each point in time, which are variables that are key inputs for the calculation of radiative transfer and precipitation processes. The Met Office Unified Model (Brown et al., 2012) makes use of the prognostic cloud Abbreviations: PC2, prognostic cloud prognostic condensate; PDF, probability distribution function. fraction and prognostic condensate scheme (PC2) in its global configuration (Wilson et al., 2008). This is used for both global forecasting and global climate simulations. The PC2 scheme is also used in the tropical configuration of the model used for regional atmospheric modelling (Bush et al., 2019).
The parameterisation schemes in the Met Office Unified Model (MetUM) are under continuous development in order to incorporate new physical understanding, address long-standing model biases or to improve the coupling and consistency between different schemes. The thermodynamic properties, cloud cover and condensate amounts calculated by PC2 for example, imply a certain thermodynamic variability, which could provide useful inputs to other parts of the moist physics. This variability information consists of two parts, representing the clear and cloudy portions of the gridbox. The manner in which the two parts are blended together at the saturation boundary has implications for the magnitude of cloud fraction changes calculated by the PC2 scheme. Additionally, if any other physics scheme wishes to make use of sub-grid thermodynamic information, it would be best if it used the same assumptions as are used in PC2, in order to avoid inconsistencies between schemes. However, evidence of undesirable aspects of the current blending method is presented. A new blending method is proposed, which removes the unwanted behaviour. Since the impact of the existing inconsistencies is currently being overwritten before they can affect the simulation, the new formulation has negligible impact on model simulations. However, the key benefit is a more sensible basis for self-consistency between schemes in the future.
Section 2 reviews the thermodynamic variability framework used by the PC2 cloud scheme as well as describing how the offline tests and MetUM simulations were set up. Examples of poor, and improved, numerical behaviour are presented in Section 3. A discussion appears in Section 4.

| Cloud scheme framework
The PC2 cloud scheme (Wilson et al., 2008) uses the sdistribution framework, developed by Mellor (1977) and Sommeria and Deardorff (1977). Some of the key details, taken from Wilson and Gregory (2003), are reproduced here. The total humidity, q T , is defined as the sum of the vapour, q v , and condensed liquid water contents, q cl . The grid-box mean liquid water content is then defined so as to prevent any super-saturation, that is, An overbar denotes a grid-box mean and q sat is the saturation specific humidity at temperature T and pressure p. Following Wilson and Gregory (2003), the local liquid water content can be rewritten as where Q c = a L q T −q sat T L ,p À Á À Á describes the mean properties of the grid-box and s = a L q T 0 −αT 0 L −βp 0 À Á represents the deviations from that mean. Following the nomenclature in Smith (1990) and Wilson et al. (2008): and β is the derivative of q sat with respect to pressure, keeping temperature constant. Assuming that G(s), the shape of the s probability distribution function (PDF) is known, then the liquid water content, q cl , and cloud fraction, C, can be calculated from it using: and Still following Wilson and Gregory (2003), changes to the liquid water content and liquid cloud fraction in PC2 are calculated using: and where As a result, provided there is partial cloud cover and condensate present, then changes to them can be calculated without knowing all the details of the s PDF, all that is required is G(−Q c ), the height of the s PDF at the boundary between the clear and cloudy air. Calculating changes to the cloud variables in this way is known as applying "homogeneous forcing," the assumption being that the temperature, humidity or pressure forcing in Equation (10) is applied homogeneously over the whole grid-box. As a result, the forcing only modifies Q c , the grid-box mean value, the shape of the s PDF is unchanged. Wilson and Gregory (2003) assume that the PDF consists of two parts described by polynomials. The righthand, cloudy, end of the PDF is given by: while the left-hand, clear, end is given by: where the saturation deficit, SD, is given by: The parameter n describes the order of the polynomial; n = 0 implies a top-hat distribution, while n = 1 implies a triangular shape. In order to calculate the change in the cloud fraction using Equation (9), one needs to determine G(−Q c ), the height of the s PDF at the saturation boundary. This done by blending the two solutions, that is, Following Wilson and Gregory (2003), this was implemented into PC2 (Wilson et al., 2008(Wilson et al., , 2008b using: with m = 0.5, such that the right-hand side PDF dominates when the cloud fraction is small and the left-hand side PDF dominates when the cloud fraction is large. Figure 1 illustrates an example situation. Given a thermodynamic state with temperature, humidity and pressure as given in Figure 1a, the moister (blue) and dryer (brown) parts of the PDF will be as shown. Note that the s PDF has a mean of zero and the saturation boundary separates the clear and cloudy parts of the PDF. Note that as in PC2, top-hat distributions (n = 0) are used here. Figure 1b shows the impact of 0.22 K of cooling, which moves the saturation boundary to the left, meaning that more of the PDF will be above saturation. Figure 1c shows how the cloud fraction change is calculated as the product of the change in the saturation boundary, ΔQ c , and the height of the PDF at that boundary. That height, G(−Q c ), is itself calculated as a blend of G 1 and G 2 , the heights of the PDF in the cloudy and clear parts.
From Figure 1 and Equations (7) and (11), it is worth noting that the right-hand side PDF has an area proportional to C, which is the product of its height, G 1 / C 2 q cl , and its width, / q cl =C . In the limit of dissipating cloud (C ! 0), the combined PDF, G(−Q c ), should tend towards the G 2 , the PDF in the clear-sky. However, the size of G 1 , will depend on how C and q cl tend to zero. Given that C and q cl can decrease independently, it is possible to have finite C with very small q cl . The very small denominator in the definition of G 1 can then make G 1 very large indeed.
Closing the G(−Q c ) equation using a blending based on the area of the two PDFs leads to these very large values of G 1 affecting the value of G(−Q c ). One way of ensuring better numerical behaviour, and ensuring that in the limit of clear sky G(−Q c ) ! G 2 , is to reduce the F I G U R E 1 An example of the thermodynamic probability distribution function (PDF) framework used in the Met Office Prognostics Cloud scheme (PC2): (a) the cloudy (moister) part of the PDF, on the right (cyan), has height G 1 = C 2 =2 q cl and area C; the clear (dryer) part, on the left (brown), has height G 2 = (1− C) 2 /2SD and area 1 − C, the mean of the PDF is zero; (b) the impact of a temperature change of −0.22 K on the new position of Q c (dashed line); (c) the change in the cloud fraction is shown by the green shading, found from the product of the width of the Q c change and the height G(−Q c ), a blend of G 1 and G 2 . Note that unlike Figure 5, the x-and y-axes all use a linear scale weight of the contribution from G 1 as clear sky conditions are approached. This can be achieved by blending the PDFs using weights based on their widths, rather than their areas.
In the other limit, as one approaches overcast conditions, G(−Q c ) should tend towards G 1 . However, if the SD is very small, G 2 can become very large. One way of mitigating this, is for the contribution of G 2 to G(−Q c ) to reduce as C ! 1.
In summary, rather than blending the PDF contributions based on their areas, it is proposed to blend them in proportion to their widths. Specifically, it is proposed to redefine the weights as:

| Idealised offline tests
Offline tests are performed to explore the parameter space and find where unrealistic values of G(−Q c ) are being produced. To generate a series of q cl and C states from which to calculate G(−Q c ) offline, we consider a grid-box at T = 5 C and p = 800 hPa. For these tests the relative humidity, RH crit , at which cloud first starts to form is assumed to be 0.8. Then, assuming a symmetric triangular PDF and considering a range of total relative humidities, RH T = q T =q sat T L , p À Á , between 0.8 and 1.2 (80-120%), based on Wilson et al. (2007), we define a series of q cl and C initial states using: and where Þis the width of the PDF. Each element of the q cl vector is then associated with each of the elements of the C vector to generate a twodimensional parameter space of plausible q cl and C pairs, from which we can calculate G 1 , G 2 and the two ways of F I G U R E 3 An example where the dry PDF is much higher than the moist one. Panels (a-h) are maps of a 4 × 4 region centred on 41 N41 W (between Newfoundland and the Azores) at a height of 1,300 m at 1200UTC (T + 12) on January 1, 2017. (a) cloud fraction, (b) log 10 liquid water content, (c) temperature increment, (d) saturation deficit, (e) original cloud fraction change, (f) new cloud fraction change, (g) log 10 height of left-hand side PDF and (h) log 10 height of righthand side PDF blending to find G(−Q c ). Assuming the simple case where there is only a temperature forcing (and no moisture or pressure changes, hence simplifying Equation (10)), Equations (9) and (10) are used to find the smallest |ΔT| that would lead to unrealistic situations where |ΔC| > 1.

| Met Office Unified Model simulations
To provide data for investigating whether the calculation of G(−Q c ) could be leading to occurrences of |Δ C| > 1 within a model simulation, a 12-hr MetUM simulation was carried out using the GA6.1 configuration (Walters et al., 2017). The initial conditions were provided by an analysis valid at 0000Z on January 1, 2017 and the resolution was N768. The global MetUM uses a regular latitude-longitude grid with 2 × N points of longitude and 1.5 × N points of latitude. Temperature, pressure, vapour and liquid water contents, and liquid cloud fraction were output from this run in order to again calculate G 1 and G 2 offline. The two ways of blending to find G(−Q c ) were then compared.
To assess the impact of changing the PDF blending approach, a series of 5-day global weather forecasts were run. For these, the resolution was N320 and 24 forecasts (with start-dates spread through-out the annual cycle) were run and evaluated against analyses and sondes using the standard weather verification process used for model development (e.g., Walters et al., 2017). To assess the impact of the new blending formulation when running the model in climate mode, a 10-year atmosphereonly N96 simulation was performed.
It is worth highlighting that within the PC2 code, Equation (9) is actually used to find the new updated cloud fraction using: The cloud fraction increment that is output as a diagnostic is then recalculated from: This means that although the value of G(−Q c ) could be implying very large |ΔC|, the code prevents this leading to unphysical values of C and the diagnostic coming out of the model reflects the increment that was actually applied, rather than the one that was initially calculated. Although, this means that unrealistic values are not allowed to affect the model run, it does mean that the undesirable behaviour cannot be seen in the diagnostics and that offline tests, such as the ones presented here, are required to find evidence of very large |ΔC| being calculated.  Figure 2b by contrast shows G(−Q c ) for the new blending formulation. The peak value is now much smaller, of the order of 1.2e3. The shape of the G(−Q c ) appears better behaved across the parameter space, tapering to lower values at low liquid water content, rather than rising sharply. Figure 2c,d shows the magnitude of temperature change that would lead to a cloud fraction change of 1. In the original formulation, there are many parts of the parameter space, where modest temperature changes, of the order of 0.1 K over the timestep, would lead to the grid-box going from partially cloudy to completely overcast or clear within one timestep. In the new method (Figure 2d), the system is more stable; temperature changes over a single timestep would need to be larger than 5 K over large parts of the parameter space for large cloud fraction swings to occur. The minimum temperature change in Figure 2d is 3.7 K, since this is quite a large temperature change to occur over one timestep, the new blending is much less likely to lead to large cloud fraction changes which would then need to be capped.

| Offline results
T A B L E 1 Values of input cloud, thermodynamic and forcing parameters used to output cloud fraction changes offline at two example locations at 1200UTC (T + 12) on January 1, 2017 using data from a N320 global simulation

| Met Office Unified Model results
Using model data taken from a height of 1,300 m above the model surface, G 1 and G 2 are calculated and the locations of the largest and smallest ratio between them are found. Maps of the 4 × 4 regions around those two points are plotted, focussing on the variables of interest. Figure 3 shows the example where the dry PDF is much taller than the moist PDF. Panels (a) and (b) show C and q cl . As a forcing term for these offline tests, we take the long-wave radiative heating. This is shown in Figure 3c. The saturation deficit (Equation (13)), is shown in panel d). The cloud fraction change calculated using the original and new methods F I G U R E 4 As Figure 3 but for a case where the moist PDF is much higher than the dry one. The maps are centred on 27 N73 W (between the Bahamas and Bermuda) at a height of 1,300 m at 1200UTC (T + 12) on January 1, 2017 are shown in Figure 3e,f, while the numerical values of the key variable at the central location of these maps are given in Table 1. The dry PDFs (brown) are shown on the left and the moist ones (cyan) on the right. Note that unlike Figure 1, here a logarithmic scale is used for both the height and width. The horizontal lines show the blended value using the original (solid red) and new (dashed green) G (−Q c ) closure formulations F I G U R E 6 Set of three log10 two-dimensional normalised histograms, showing the distribution of data within the liquid water content, saturation deficit and cloud fraction parameter space. Additionally, the location of grid-boxes where setting G(−Q c ) = G 1 (−Q c ), or G (−Q c ) = G 2 (−Q c ), would lead to a cloud fraction change of more than 2 in magnitude, are shown by circles, or crosses, respectively F I G U R E 7 Equitable threat score, as a function of lead time, for total cloud amount greater than 0.3 assessed against surface observations; (a, b) northern hemisphere, (c, d) the tropics, (e, f) southern hemisphere; (a, c, e) for JJA, (b, d, f) for DJF F I G U R E 8 Profiles of mean temperature error at T + 120 hr compared to sondes; (a, b) northern hemisphere, (c, d) the tropics, (e, f) southern hemisphere; (a, c, e) for JJA, (b, d, f) for DJF occurrences of values less than −1. For the point of interest in the centre of these maps, despite a modest warming of 0.03 K, an unphysically large cloud fraction change of −15 is being calculated. This is despite the C, q cl , SD and ΔT all being reasonable and not very noisy. Maps of G 2 and G 1 are shown in panels (g) and (h). A graphical representation of the two PDFs is shown in Figure 5a,b, along with the values of G(−Q c ) from the two blending methods. Note that unlike in Figure 1 where the PDFs were plotted in linear units, here, due to the vastly different height of the two PDFs, their heights are shown in log 10 terms. Similarly, due to the vastly different widths, these are also shown as log 10 of their absolute width, and a reversed x-axis is used to depict the left-hand PDF in Figure 5a,c. From Table 1, we see that in this case, the SD is very small, despite the cloud cover not being that close to 1. This makes G 2 several orders of magnitude larger than G 1 , and with the original weighting given by Equation (15) this gives a final G(−Q c ) shown by the solid red horizontal line in Figure 5a,b. The new weighting developed here (Equation (16)), pulls the blended G(−Q c ) away from the larger PDF and gives a blend shown by the dashed green line in Figure 5a,b. The cloud fraction change associated with ΔT = 0.03 is now a more sensible −0.007, instead of −15 (Figure 3f and Table 1).
The second example, Figure 4, has the opposite scenario in terms of the PDFs, with the moist PDF much higher than the dry (Figure 5c,d). This time a plausible ΔT = −0.15 is leading to the cloud fraction change by +11. With the new blending, G(−Q c ) is greatly reduced and this becomes a more reasonable +0.01 change in cloud fraction (Table 1).
Having looked at the example of those two grid-points in detail, Figure 6 considers the whole of the layer at 1300 m altitude. Considering the variables of q cl , C and SD, two-dimensional histograms are calculated. The normalised frequency of points in these parameters spaces is shown, in log 10 terms, by the shading on the 3 faces of Figure 6. For the q cl and C parameter space, we see that the main body of points follows a pattern where larger water contents are generally associated with larger cloud fractions. However, and this was one of the original aims of PC2 (Wilson et al., 2008), a given water content can be associated with almost any cloud fraction and similarly a given cloud fraction can have a range of water contents spanning more than an order of magnitude. The locations where the magnitude of the G 1 or G 2 terms (if taken to represent the value of G(−Q c ) on their own) when combined with the long-wave heating rate would lead to a cloud fraction change of more than 2 in magnitude, are shown by black circles (for G 1 ) or black crosses (for G 2 ) in Figure 6.
Large G 1 occurs over a range of cloud fractions and, as would be expected from Equation (11), are always associated with low water contents that are within the extreme low end of the q cl population. Similarly, G 2 is associated with very low values of saturation deficit, at the extreme end of that population. Very large values of G 1 and G 2 can occur over the full range of cloud fraction values.
Physically, these two extremes of PDF behaviour can be explained as follows. If one considers the scenario of stratocumulus, where q cl can be depleted as it is converted to rain and precipitates out, while C remains steady, then the width of the right-hand side PDF will reduce (in proportion to the reduction in q cl ) while G 1 , the height of the PDF, will increase to ensure that the product of width and height leads to a constant cloud fraction. Similarly, for the left-hand side PDF, if the gridbox mean vapour amount approaches saturation while the sky is not overcast, then this means that the width of that PDF is very small, while the height must increase to maintain a product that is equal to (1 − C).
Concerning how often unphysical increments are occurring and being capped, we note that, in these N768 simulations, each model level consists of around 1.7 million points. For the instantaneous model field at an altitude of 1,300 m considered here, around 25% of grid boxes (430,000) had liquid cloud in them. Of these, around 250 grid boxes have unphysical cloud fraction increments whose magnitude is larger than 1. Although around 0.05% is a small fraction of the cloudy points, it does suggest that over the course of a day of simulation, unphysical increments are happening several thousand times.
In terms of the impact of implementing the new blending method, since the change to the code directly affects cloud cover we consider Figure 7. This shows the equitable threat score for cloud cover greater than 0.3, as a function of lead time for a series of weather forecasts. The error bars encompass a range of one standard deviation around the mean. Irrespective of the region or time of year, there are no statistically significant changes in the characteristics of the cloud cover error.
T A B L E 2 Globally averaged, 10-year mean, short-wave (SW) and long-wave (LW) cloud radiative effect (CRE) in the CERES EBAF observations and in the GA7 control (CTL) and in the experiment (EXPT) with the new blending method As a measure of indirect impact on the manner in which the cloud fraction is updated, Figure 8 shows the changes in the temperature error, when assessed against radio-sondes for a lead-time of 120 hr. Again, there are no significant changes in the model performance.
With regards to implementation in a climate simulation, Table 2 summarises the global-mean cloud radiative effects from the 10-year atmosphere-only climate simulation. There is a small reduction in the biases and RMSE, but essentially it is possible to fix the inconsistency without deteriorating model performance.

| DISCUSSION
It has been shown, using offline tests and analysing some actual MetUM data, that the current method of blending the two PDFs in PC2 in order to calculate homogeneous forcing increments for cloud fraction can lead to unphysically large cloud fraction changes.
The NWP case studies and AMIP simulation show that the new way of blending the PDF to find G(−Q c ) has no significant impact on the model evolution. Although this lack of beneficial impact could be disappointing, it is understandable since the PC2 code had been written including a cap to ensure that the cloud fraction could never go outside of the range of zero to one. As a result, the unrealistically large cloud fraction changes were never allowed to push the model into an unphysical state. Interestingly, as the cap is applied immediately prior to the diagnostics being calculated and output, it was hard to know that there was an issue, until detailed analysis such as the one presented here was carried out. This analysis does highlight the importance of checking the numerical behaviour of model parameterisation schemes, not just by looking at the output diagnostics, but also by looking at how they behave offline using plausible test data.
In terms of future work, the shape of the PDF, and the value of G(−Q c ), could be validated by comparing the PDF inferred by PC2 in a coarse-resolution model to some high resolution simulations done for the same time and place. This could be done using 100-m simulations within the MetUM (Hughes et al., 2015;Hanley et al., 2016) or one could make use of the large-eddy simulations performed for a whole year by Schalkwijk et al. (2015), which have the added advantage of being located over a well instrumented observational site.
To conclude, the neutral impact on model performance is reassuring, as it allows the new formulation to be implemented easily, without needing to retune the model. The new formulation then gives the benefit that it will now be possible for other physics schemes to make use of sensible thermodynamic PDF information, and for that information to be consistent with what the PC2 scheme is doing. For example, this could include using thermodynamic PDF information when deciding whether or not a convective parameterisation should be triggered. It will also be relevant when calculating how subsidence warming associated with convection is affecting the large-scale state and any pre-existing clouds within it.