Demystifying the complex nature of stratocumulus clouds with machine learning

Why do we need to model stratocumulus clouds? Clouds reflect the Sun’s energy and insulate the Earth by absorbing longwave radiation, making them a key part of our planet’s energy budget. Stratocumulus clouds are shallow and usually located in the lower 2km of the atmosphere (Figure 1). They cover around a fifth of the Earth’s surface, with extensive near-permanent cloud fields located to the west of all major continents (Figures 2 and 3). This large coverage over the darker land or ocean surface combined with bright, reflective cloud tops acts to increase the reflectivity of the planet (low cloud cover alone increases the Earth’s albedo by ~25%, from about 0.24 to 0.30 (Hartmann et al., 1992)), whilst emitting longwave radiation with approximately the same energy as the surface. The result is a net cooling effect in most regions where these clouds are present (Wood, 2012). Stratocumulus clouds are impacted greatly by changes in aerosol concentrations. Aerosols are gas, liquid or solid particles that are suspended in the air. Aerosols can come from natural sources, such as ocean spray or deserts, or anthropogenic sources, like smoke from vegetation burning and industrial sites. Aerosols directly affect the radiation budget through absorption and scattering of solar energy, but they also indirectly affect it through their interactions with clouds. They are key to cloud formation because they provide a surface onto which water vapour in the air can condense, at a relative humidity very close to 100%, to produce a cloud droplet. Without aerosol particles a much higher relative humidity (around 300%) would be required. As such, all cloud droplets contain aerosol particles and the properties of the cloud, such as brightness and cloud cover, are altered when more aerosol enters or leaves the cloud. These processes are known as aerosol– cloud interactions (Bellouin et al., 2020). Considering the importance of stratocumulus clouds and their cooling effect, Demystifying the complex nature of stratocumulus clouds with machine learning


Why do we need to model stratocumulus clouds?
Clouds reflect the Sun's energy and insulate the Earth by absorbing longwave radiation, making them a key part of our planet's energy budget. Stratocumulus clouds are shallow and usually located in the lower 2km of the atmosphere (Figure 1). They cover around a fifth of the Earth's surface, with extensive near-permanent cloud fields located to the west of all major continents (Figures 2 and 3). This large coverage over the darker land or ocean surface combined with bright, reflective cloud tops acts to increase the reflectivity of the planet (low cloud cover alone increases the Earth's albedo by ~25%, from about 0.24 to 0.30 (Hartmann et al., 1992)), whilst emitting longwave radiation with approximately the same energy as the surface. The result is a net cooling effect in most regions where these clouds are present (Wood, 2012).
Stratocumulus clouds are impacted greatly by changes in aerosol concentrations. Aerosols are gas, liquid or solid particles that are suspended in the air. Aerosols can come from natural sources, such as ocean spray or deserts, or anthropogenic sources, like smoke from vegetation burning and industrial sites. Aerosols directly affect the radiation budget through absorption and scattering of solar energy, but they also indirectly affect it through their interactions with clouds. They are key to cloud formation because they provide a surface onto which water vapour in the air can condense, at a relative humidity very close to 100%, to produce a cloud droplet. Without aerosol particles a much higher relative humidity (around 300%) would be required. As such, all cloud droplets contain aerosol particles and the properties of the cloud, such as brightness and cloud cover, are altered when more aerosol enters or leaves the cloud. These processes are known as aerosolcloud interactions (Bellouin et al., 2020).
Considering the importance of stratocumulus clouds and their cooling effect, Demystifying the complex nature of stratocumulus clouds with machine learning  research is needed to reduce the uncertainty in aerosol-cloud interactions because it dominates the uncertainty in estimates of the total radiative forcing. In addition to aerosol-cloud interactions, stratocumulus clouds can be altered by changes in the climate, and this could have consequences that feed back into the atmosphere. For example, the increase in CO 2 emissions warms the atmosphere and consequently the oceans. Current models suggest that warmer surface temperatures will reduce the amount of low, shallow clouds, which will decrease the planet's reflectivity and allow more solar energy to be absorbed (Ceppi et al., 2017). This feedback would act to amplify the warming from the original CO 2 emissions. Precise predictions of our future climate rely on our understanding of how stratocumulus clouds will respond to changes in climate and what climate feedbacks might occur.

Why are clouds difficult to model?
Knowledge of aerosol-cloud interactions in shallow clouds has increased rapidly over the last few decades due to increased observational campaigns and satellite data. Recent observations provide clear evidence for some effects, such as the Twomey effect (Twomey, 1977;Malavelle et al., 2017). For a cloud with an unchanging volume of liquid water, an increase in aerosol particles will distribute the water volume between more droplets, creating smaller and more numerous cloud droplets. The adjusted cloud is brighter and reflects more sunlight. Also, the lifetime effect may occur in clouds that are likely to precipitate. Here, a higher aerosol concentration will suppress precipitation because smaller droplets collide less efficiently. Retention of moisture grows the cloud and prolongs the cloud lifetime (Albrecht, 1989;Christensen et al., 2020). Observations are critical because they provide evidence (or lack thereof ) for real-world effects and they are a benchmark for model comparison.
Other aerosol-cloud interactions are less clearly observed because they are dependent on meteorology, radiation and the properties of the surface below. Following the Twomey effect, the observed outcome could be cloud growth or cloud dissipation. Figure 4 shows the key stratocumulus processes, which are balanced for a stable cloud layer. Heating and cooling processes contribute to the movement of air (turbulence) and evaporation or condensation of cloud droplets. Turbulence at cloud top causes entrainment, the folding of dry air from above the cloud into the moist cloud air, which evaporates cloud droplets. If processes that increase turbulent mixing and moistening in the cloud dominate, the cloud grows. If processes that remove aerosol (drizzle) or increase evaporation (entrainment) dominate, the cloud dissipates. There are many aerosol-cloud interactions and each one influences the balance of these processes. Often this results in multiple processes compensating each other and a stable cloud is maintained (Stevens and Feingold, 2009). But stratocumulus are also susceptible to sudden changes in behaviour when conditions are met that surpass certain thresholds.
A major challenge in modelling clouds within the climate system is the range of scales that need to be covered with finite computational resources (Mülmenstädt and Feingold, 2018). Global climate models simulate large atmospheric circulations over thousands of kilometres with grid boxes typically covering a few hundred kilometres, much larger than cloud-scale eddies that occur on the scale of tens of metres. Instead, the processes that cannot be explicitly resolved are represented by approximations called 'parameterisations' , but there are often multiple ways to achieve them. The discrepancies between models that use different parameterisations create a large range of model predictions for the future climate. Detailed studies of cloud processes, along with observations, are required to improve parametrisations to find a convergence point that reduces the uncertainty in climate predictions.
For every physical process described in a model, there are inputs, or 'parameters' , for which we do not know the exact value and it could be any number within some range, making it uncertain (Carslaw et al., 2018). The uncertain parameters are constants, such as droplet terminal velocity, or atmospheric conditions, such as wind speed. If a process has two uncertain parameters, P1 and P2, each with a possible range of 0-10, there is a two-dimensional parameter 'uncertainty space' (henceforth, parameter space) of possible combinations of values. To simulate cloud behaviour at each integer point in this parameter space would take 100 simulations. The 'curse of dimensionality' means that an increase in the number of dimensions (here, the number of uncertain parameters) will require exponentially more simulations to get an adequate coverage of the parameter space. Considering the many uncertain parameters in modelling clouds, we need a feasible method for exploring parameter spaces.

Statistical emulation: a vital tool for modelling clouds
One tool for exploring parameter spaces comes from machine learning using Bayesian statistics. In classical statistics, probability is viewed as the frequency of an event happening given infinite repetitions of an action. The uncertainty here comes from the inability to make exact predictions when there is some randomness. For events that are not repeatable, uncertainty comes from the gaps in our knowledge about the event (O'Hagan, 2006). For example, suppose we want to throw some different-sized balls off the top of a building and find the average fall time. We make an initial guess and after throwing one ball, we update it based on what we find. Each time we throw one, we update our guess and narrow the uncertainty. Bayesian statistics allows for this type of uncertainty because probability can be more widely classified as the confidence that an event will occur based on prior experience. By following a Bayesian approach, where we update our prior belief based on new information, an informed prediction on the next outcome can be achieved.
Using Bayesian statistics in computer modelling greatly reduces the number of model simulations required to explore parameter spaces (O'Hagan, 2006). Any system of processes that we want to study can be represented by an unknown function that describes the transformation from the input parameters to a corresponding output of interest. Using a method called statistical emulation, the function can be approximated. It starts with a prior specification based on any initial knowledge of the system. Then, a set of simulations can inform us about the model's behaviour over parameter space by simultaneously perturbing all parameters across the ranges of their possible values. This set of 'training data' is used to update the prior into a better approximation of the unknown function. The updated approximation, known as the posterior distribution, is an emulator for that specific output of interest and it can be used to make predictions at any new combination of parameter values within the parameter space.
Statistical emulation has been widely used to approximate complex models. A common motivation is to identify which uncertain parameters contribute most to the model output uncertainty, called sensitivity analysis. This has been used in studies on aerosol properties from global models (Lee et al., 2011;Carslaw et al., 2013;Regayre et al., 2014). Emulation is also used alongside 'history matching' to rule out values in the parameter ranges that contradict real world observations (Williamson et al., 2013;Andrianakis et al., 2017). Recently, emulation of cloud models was used to study the influence of environmental and aerosol conditions on cloud properties (Wellmann et al., 2018;Glassmeier et al., 2019;Hoffmann et al., 2020;Park et al., 2020). The ability to study how changing multiple parameters at once affects the output across the full range of uncertain values is particularly useful in modelling the many interacting cloud processes.
The following case study is an example of how statistical emulation can be applied to cloud modelling. We consider just two uncertain parameters for the ease of visualising the results, but this method can be scaled up to tens of parameters, including aerosol concentrations and processing parameters to study aerosol-cloud interactions. The remainder of this article will discuss simulation set-up, results from the cloud model, and the emulator results.

Simulation set-up
Large-eddy simulation models are run at a high enough resolution to resolve the largescale dynamics of cloud formation and development. Typical resolutions in largeeddy simulations of clouds are in the range of tens to hundreds of metres, horizontally, and 5-20m vertically. Many processes are still parameterised, such as droplet interactions, small-scale turbulent mixing, and heating/ cooling at cloud top ( Figure 4). However, these models, informed by observations, are important tools for studying cloud processes in detail. Here, we used a large-eddy simulation model from the UK Met Office and the Natural Environment Research Council, called the MONC (Met Office/NERC Cloud) model (Brown et al., 2020). MONC is used to study dynamics and, when coupled with the Cloud AeroSol Interactive Microphysics scheme, aerosol-cloud interactions (Böing et al., 2019;Poku et al., 2019). The simulation duration was 8 hours and covered a horizontal domain size of 8.75km 2 with a resolution of 35m and a variable vertical resolution between 7.5-10m up to 1.5km.
The base simulation was representative of the observations from a research campaign in July 2001 off the Californian coast, called the Second Dynamics and Chemistry of Marine Stratocumulus field study (DYCOMS-II) (Stevens et al., 2003). Research flight RF01 recorded a uniform stratocumulus field with no drizzle, and a well-mixed boundary layer capped by a temperature inversion. The boundary layer is the lower portion of the troposphere that is most affected by surface winds and heating, where cloud-forming convection is initiated ( Figure 4). When it is well-mixed, there is an even overturning of air throughout the layer which keeps the defining properties mostly constant throughout. The top of the boundary layer is often capped by a temperature inversion, above which the properties of the free troposphere are different. The conditions there have a strong influence, via entrainment, over whether stratocumulus cloud will be maintained and grow or whether it will transition into cumulus clouds or disappear entirely.
To study the transition between cloudgrowing and cloud-dissipating conditions, we perturbed two free-tropospheric properties from the base simulation. Both perturbed parameters are related to the change in conditions at the inversion. The first parameter was the size of the increase in potential temperature, ∆θ, and the second parameter was the size of the decrease in total water mass-mixing ratio, ∆q t : That is, changing how much warmer and drier the free-tropospheric air was than the cloudy air into which it is entrained. Figure 5 shows the profiles for these parameters through the boundary layer in dark green for the base simulation, and the inversion at the free troposphere is indicated by the jump at 840m. Uncertainty ranges for the parameters were taken from a theoretical study that derives the relationship between these two parameters and cloud thinning or thickening (Van Der Dussen et al., 2014). In this study, we consider a measure of the liquid water content integrated through the boundary layer (an indicator of cloud growth) called the liquid water path (LWP).
A Latin hypercube method was used to choose which combinations of values in parameter space to simulate with the model. The method has been shown to be optimal for statistical emulation because it designs a set of simulations that covers the twodimensional parameter space well without overlapping in each one-dimensional range (Loeppky et al., 2009). The final design for our training data can be seen, in blue, to be well-spaced in two dimensions (Figure 6(a)) as well as in one dimension (best seen in Figure 5).
An extra set of simulations are run as validation data, shown in pink in Figure 6(a). These data are not used to train the emulator but are used instead to check the emulator is an accurate approximation of the model. The total number of simulations for this case study was 33, consisting of the base, four tests to check the extremes of parameter space, 20 training and eight validation simulations. This case study will be described in detail in a subsequent study.

Model results
Figure 6(b) shows the training data LWP, averaged over the last two hours of the 8-hour simulation. These results show a clear pattern across parameter space. At ∆q t = 0 and ∆θ = 0, the air above the cloud would be the same as the cloudy air, so every point in parameter space has warmer, drier abovecloud air than in the cloud. The coolest, driest conditions are the bottom-left corner of parameter space and the warmest, wettest conditions are the top-right corner. LWP increases with increasing ∆q t , from a minimum of 6gm −2 to a maximum of 176gm −2 .
LWP also generally increases with ∆θ, but this is much clearer when ∆q t is more negative than when near to 0. The behaviour in parameter space shows very thin clouds in region A, but with the increase in LWP into region B, thick, drizzling stratocumulus were simulated. The physical change in cloud behaviour across the parameter space can be interpreted as follows.
All the simulations in region B were drizzling from the beginning of the simulation and this process dominated. As the drizzle fell from the saturated cloudy air into the sub-saturated air below, much of it evaporated and cooled the sub-cloud air. The drizzle also removed water from the cloud-top region which prevented the cloud from growing upwards with convection. The cloud tops in most of these simulations either stayed at the same height during the simulation or sank a little in the area with higher ∆θ (stronger temperature inversion). The cooler sub-cloud air encouraged more condensation of cloud droplets in this area of parameter space which created convection and cloud growth in the sub-cloud layer. Hence, these clouds had a high LWP because they grew into the sub-cloud layer and thickened.
There was no drizzle throughout the simulations in region A. The free-tropospheric conditions were just slightly warmer but very dry compared with the cloud layer. Entrainment of this air into the cloud strongly evaporated cloud droplets and cooled the surrounding cloud top more. This cooling created more turbulence and therefore more entrainment and evaporation, forming a feedback loop (Mellado, 2017). Energetic mixing penetrated the weak temperature inversion cooling the air just above the cloud, so the inversion height rose, and the boundary layer deepened. Additionally, evaporative cooling in the cloud layer stopped the boundary layer being wellmixed, which cut the cloud off from the ocean's water source. This, combined with rapid evaporation from above, prevented cloud growth so these clouds were high and thin, resulting in a lower LWP.

Emulation
We used our training data to create an emulator that approximates the unknown function between the two parameters, ∆θ and ∆q t , and the output of interest, LWP. We used the validation data set to evaluate the emulator. Figure 7 shows the model's output at the eight validation points against predictions made using the emulator for the same points in parameter space. The black line is the line of equality where y(x) = x and the error bars show the 95% confidence bounds on the emulator predictions. All the markers are close to the black line and this line falls well within the confidence bounds of the predictions, so we infer that this emulator is a good representation of the relationship between the two parameters and the output. Hence, we can confidently use predictions from this emulator to study the LWP response to changes in ∆θ and ∆q t .
Here, we are specifically using Gaussian process statistical emulation. It has one key advantage over other machine learning techniques, in that each predicted point has a calculated Gaussian error that is associated with using the emulator instead of the complex model itself. So, each predicted value for LWP has a 95% confidence interval. Close to the training data this error is zero because we know the model's output at that point, but predictions made at a point far from the training data have a larger associated error.
We interpolated over the whole parameter space using the emulator's estimated function to make predictions for average LWP at 2500 new combinations of ∆θ and ∆q t values, within our chosen ranges. This is shown in Figure 6(c) as a response surface, which visually shows the relationship between inputs and output. The response surface shows that the transition between cloud-growing and cloud-dissipating regions of parameter space is a smooth, but quite steep, slope. One of the key assumptions when emulating is that the model output is smooth, so the potential for sudden changes in cloud behaviour can cause problems. In this case, the output is smooth enough that a standard emulation method works.
The response surface confirms what we expected from the training data in Figure 6(b) but gives a lot more information to use in further analysis. The slope through the centre of parameter space marks a change in cloud behaviour from drizzling, cloud growth into the sub-cloud layer to deepening of the boundary layer with rising, evaporating clouds. The surface has some bumpiness on it due to the natural variability of clouds, which is reflected in the model. The bumps can be characterised and smoothed in a physically justifiable way within the statistical model. This process will be studied in detail in a subsequent paper.

Concluding remarks
In this case study, we have demonstrated the use of statistical emulation to approximate the relationship between a two-dimensional parameter space and an output of interest for an idealised cloud simulation. We used a high-resolution cloud model to create a base simulation, from which we perturbed two parameters across reasonable ranges. This training data was used to create an emulator that accurately represented the relationship of interest from the model and allowed us to make predictions at thousands of new parameter combinations. The resulting response surface informs us about the cloud behaviour in different regions of the parameter space and the transition between regions. This response surface can be used for further analysis without running more simulations with the cloud model.
Each of the 33 simulations used in this case study took between 7 and 10 hours of realworld time to run on a high-performance computer. So, in total, this small number of simulations took more than 250 hours of computing time. Predicting the output at 2500 new points in parameter space took less than a second on a standard laptop, and the associated error is quantified. This is a small example of how emulation can be used. However, it is also a demonstration that machine learning tools are invaluable for studying the behaviour of these cloud systems. Other machine learning techniques, such as neural networks, are excellent for applications like classification of cloud structures since we have large quantities of satellite data (Denby, 2020;Juliano and Lebo, 2020). But in modelling studies where the quantity of data may be limited by computational resource, a method like statistical emulation is more appropriate.
As mentioned, in statistical emulation we assume the model's output is smooth, which could be a problem when studying cloud systems. Smooth means that a small change to an input will make only a small difference to the output. Cloud behaviour, and therefore cloud model output, is not always smooth for two reasons. Firstly, clouds are slightly chaotic in that a small change at the beginning of cloud development can have knock-on effects that result in a large change later. Secondly, when a certain threshold is met a cloud may exhibit a sudden change in behaviour. If the change is too sharp, we cannot accurately  emulate the relationship with this standard emulation method. Instead we would have to adapt the method, such as by splitting up parameter space into regions and emulating them separately or by adapting the prior specification (Mohammadi et al., 2019;Pope et al., 2021).
Currently, statistical emulation is used to explore much larger parameter spaces with up to around 50 dimensions of interest and, in some cases, multiple outputs of interest. Further work in this project will include perturbations to aerosol concentrations and more atmospheric parameters, such as wind speeds and inversion height, to study their influence on various cloud properties. We can reduce the uncertainty in aerosolcloud interactions, and hence uncertainty in estimations of radiative forcing, through idealised cloud studies like these, which allow exploration of all the interacting processes together. By combining the perturbations of multiple parameters, we get an overview of how any output of interest varies when multiple things change at once, which is representative of a real-world environment. Additionally, cloud processes in global climate models (either parameterised or explicitly calculated) contain many uncertain parameters. By studying these cloud processes and aerosol-cloud interactions with statistical emulation and high-resolution cloud models, we can learn which parameters are most influential, constrain uncertain ranges, and identify new relationships that are difficult to derive from observations.

Further reading
For more reading on the motivation for using statistical emulation in climate and cloud models, the opinion piece cited earlier in this article gives a detailed, yet accessible, explanation (Carslaw et al., 2018). Similarly, Watson-Parris (2020) gives an overview of the current state of climate model emulation and a full description of the different types of uncertainties involved in climate modelling. A more in-depth justification is given along with some pre-emulation model results in the proposed framework of Feingold et al., (2016). A general tutorial on the mathematics and implemention of Gaussian process emulation in any field is presented in O'Hagan (2006). Two studies comparing machine learning regression techniques found Gaussian process regression to give the best results (Salter and Williamson, 2016;Schlosser et al., 2020). Several studies use various machine learning techniques to train parameterisations, either to speed up an existing one or to approximate a relationship that is unknown (O'Gorman and Dwyer, 2018;Seifert and Rasp, 2020;Chiu et al., 2021).